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Preface 


The  purpose  of  this  investigation  is  to  determine  the  energy 
associated  with  the  electromagnetic  spectrum  between  100  kilohertz  and 
100  megahertz  from  high  altitude  nuclear  EMP.  The  computer  code 
developed  to  generate  simulated  EMP  fields  (EMPFRE)  is  designed  to  be 
flexible.  The  gamma  ray  pulse  can  be  represented  as  a  collection  of 
functions  or  can  be  adapted  to  read  and  interpolate  between  a  set  of 
data  points.  Other  physical  parameters  such  as  Compton  electron 
characteristics  are  easily  identified  and  changed.  A  parametric  study 
is  conducted  to  determine  cause  and  effect  relationships  between  burst 
parameters  and  the  frequency  characteristics  of  the  EMP. 
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work.  The  insight  I  gained  from  his  advice  is  appreciated. 
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Abstract 


A  program  is  developed  to  model  the  electromagnetic  pulse  from  a 
high  altitude  nuclear  detonation.  A  Runge-Kutta  numerical  technique  is 
used  to  solve  for  the  electric  fields.  A  continuous  Fourier  Transform 
of  the  EMP  is  used  to  determine  the  frequency  profile  of  the  EMP. 
Parametric  studies  are  performed  to  determine  cause  and  effect 
relationships  between  burst  parameters  and  the  EMP  frequency  profile 
from  100  KHz  to  100  MHz.  Burst  parameters  studied  are:  gamma  pulse 
time  history,  gamma  ray  energies  from  1  MeV  to  10  MeV,  gamma  ray  yield, 
height  of  burst  from  75  Km  to  200  Km  and  intersection  angle  of  the  slant 
range  with  the  geomagnetic  field  from  90  degrees  to  30  degrees. 


I.  INTRODUCTION 


BACKGROUND 

EMP  (electromagnetic  pulse)  is  a  prompt  effect  of  a  nuclear 
explosion,  consisting  of  a  broad  spectrum  of  radio  waves.  These  radio 
waves  can  induce  electric  currents  into  unprotected  systems.  The 
currents  cause  resistive  heating  that  can  melt  electrical  contacts  or 
damage  electrical  components.  Estimates  of  the  shape  and  magnitude  of 
the  EMP  are  based  on  sparse  experimental  data  and  computer  simulations 
of  the  physical  processes  which  generate  the  EMP. 

Studies  conducted  in  the  late  1950's  concluded  that  the  EMP  from  a 
high  altitude  nuclear  explosion  would  be  confined  to  the  upper 
atmosphere  and  would  not  be  a  threat  to  syster;s  near  the  ground; 
therefore,  little  attention  was  paid  to  the  subject.  EMP  resurfaced  as 
a  problem  in  1962  when  a  1.9  megaton  bomb  detonated  248  miles  above 
Johnson  Island  was  believed  responsible  for  the  failure  of  street  lights 
in  Hawaii  (800  miles  away)  (Ref  3)*  A  subsequent  experiment  to 
characterize  the  EMP  was  inconclusive  because  the  magnitude  was  still 
underestimated  and  the  equipment  used  was  saturated  by  the  EMP.  In  late 
1963,  scientists  at  Los  Alamos  Corporation  and  the  RAND  Corporation  were 
the  first  to  explain  the  physical  processes  which  generated  the  EMP. 
The  limited  test  ban  treaty  precluded  the  experiments  that  could  be  used 
to  verify  the  theory  of  EMP  generation  (Ref  3:1009). 

"Defense  strategists  today  assume  that  a  single  Soviet  warhead 
detonated  200  miles  above  Nebraska  would  knock  out  unprotecte.i 
communications  equipment  all  across  the  United  States"  (Ref  4:1248). 
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The  projected  damage  to  military  systems  are  based  on  two  assumptions. 


First,  the  EMP  will  have  peak  fields  on  the  order  of  50,000  volts  per 
meter.  Second,  the  high  frequency  end  of  the  electromagnetic  spectrum 
has  sufficient  energy  to  induce  large  electrical  currents  into  exposed 
lengths  of  wire.  The  amount  of  energy  in  the  high  frequency  end  of  the 
electromagnetic  spectrum  depends  on  the  peak  electric  field  and  shape 
parameters  such  as  the  rise  time  (time  required  to  reach  the  peak  of  the 
pulse).  The  damage  potential  is  sensitive  to  the  EMP  spectral 
characteristics . 

SCOPE 

This  investigation  determines  the  energy  associated  with  the 
frequency  components  of  the  EMP  from  a  high  altitude  nuclear  detonation. 
This  goal  is  achieved  in  two  stages.  First,  a  computer  code  (EMPl'PE'  is 
developed  to  model  the  processes  that  create  the  EMp.  Second,  a 
continuous  Fourier  transform  is  found  for  the  simulated  EMP  and 
Parseval’s  theorem  is  used  to  calculate  the  energy  associated  with  the 
frequencies  between  100  kilohertz  and  100  megahertz. 

ASSUMPTIONS  AND  LIMITATIONS 

EMPFdE  calculates  the  electrical  field  strength  at  discrete  values 
of  time  up  to  ^  shakes  (1  shake=1E-.9  seconds).  The  late  time  electric 
fields  will  not  have  significant  impact  on  the  frequency  range  of 
interest.  The  energy  distributions  are  based  on  this  time  range. 

EMPFRE  relies  on  approximations  developed  by  Karzas  and  Latter  (Ref 
9;  for  the  currents  and  secondary  electron  density.  Discrete  numerical 
methods  are  used  extensively  including  a  fourth  order  Runge-Kutta 
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technique  to  solve  for  the  electric  fields.  Discussions  of  the  impacts 
from  these  approximations  are  presented  in  Chapter  IV. 

OVERVIEW 

The  background  theory  for  high  altitude  EMP  is  presented  in  Chapter 
II.  The  methods  used  to  model  the  processes  that  generate  EMP  are 
presented  in  Chapter  III.  Sample  results  are  presented  in  Chapter  IV  to 
characterize  the  errors  due  to  program  parameters.  Parametric  studies 
relate  burst  parameters  to  frequency  characteristics  in  Chapter  V. 


Conclusions  and  recommendations  are  made  in  Chapter  VI. 


II.  THEORY 


This  chapter  is  a  distillation  of  theoretical  work  done  by  W.  J. 
Karzas  and  R.  Latter  (Ref  9),  J.  H.  Erkkila  (Ref  5,6,7)  and  C.  L. 
Longmire  (Ref  12).  Information  from  these  sources  form  the  theoretical 
foundation  needed  to  understand  the  processes  that  lead  to  the 
generation  of  high  altitude,  nuclear  EMP.  The  primary  source  of  EMP 
from  high  altitude  detonation  is  the  result  of  the  gamma  rays.  Gamma 
rays  account  for  0.1  to  0.5?  of  the  total  yield.  Gamma  rays  entering 
the  atmosphere  produce  Compton  electron  currents.  Maxwell's  equations 
predict  the  electromagnetic  fields  resulting  from  these  currents.  The 
electromagnetic  fields  that  are  transverse  to  the  gamma  ray  direction 
propagate  from  the  absorption  region  to  the  ground.  processes  relating 
to  this  phenomenon,  are  described  and  references  are  provided  for  further 
detail. 

GAMMA  RAY  DEPOSITION 

The  number  of  gamma  rays  emitted  from  a  nuclear  detonation  as  a 
function  of  time  is  written  as: 

S^(t)  =  I  f(t)  (2.1) 

where  Y  -total  energy  of  gamma  ray  yield 

E  -mean  energy  of  the  gamma  rays  emitted 
f(t)  -time  dependence  of  the  pulse 

The  time  dependence  of  the  pulse  is  normalized: 

oo 

/f(  t)(l  .  ::  1  (2.2) 

•OD 


The  prompt  gamma  ray  emission  rate  is  proportional  to  the  neutron 
population  in  the  device.  The  gamma  ray  pulse  from  a  fission  device 
wi  1  1  rise  exponential  ly  with  the  neutron  population.  Prompt  gamma  rays 
are  generated  by  inelastic  collisions  between  neutrons  and  the  device 
materials  and  by  fission  reactions. 

The  density  of  Compton  recoil  electrons  produced  per  gamma  ray 
emitted  from  the  burst  as  a  function  of  position  is  written  as 


This  expression  assumes  the  gamma  ray  is  absorbed  by  its  first 
interaction.  Multiple  scatters  are  not  accounted  for.  The  scattered 
gamma  ray  wi  1  1  lag  behind  the  pulse  due  to  the  deflection. 

Gamma  ray  deposition  profiles  are  a  function  of  height  of  burst 
(HOB)  and  gamma  ray  energy.  Figures  2.1,  2.2,  2.3  and  2.H  are  contour 
plots  of  Compton  electron  generation  density  which  are  used  to 
demonstrate  the  dependence  on  HOB  for  1.5  MeV  gamma  rays.  Figures  2.5, 
2.fi  and  2.7  are  plots  of  Compton  electron  generation  density  along 
discrete  slant  ranges  which  are  used  to  demonstrate  the  dependence  on 
gamma  ray  energy  for  an  exponential  atmosphere.  The  impact  of  these 
parameters  on  the  EMP  is  a  topic  of  the  parametric  studies  in  Chapter  V. 
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COMPTON  ELECTRON  DYNAMICS 

The  Compton  interaction  can  be  visualized  as  a  billiard  ball 
collision  between  a  gamma  ray  and  a  molecular  electron.  Conservation  of 
energy  and  momentum  is  used  to  derive  the  following  relationships. 

Kinetic  energy  of  the  Compton  recoil  electron 


(  ')  /  ''j 

-incident  gamma  ray  energy 
=  a(l  -  COSO) 

=  Ey/^o 

-rest  mass  energy  of  tne  electron  (0.51 IMeV) 

-photon  deflection  angle 

Compton  recoil  electron  and  photon  deflection  angles 

C0T((J)  =  (Ua)TAN(£)  (2.5^ 

2 

where  $  -Compton  recoil  electron  deflection  angle 

The  Kline-Nishina  cross  sections  for  Compton  scattering  describe  the 
scattering  probability  verses  angle  (Ref  5).  The  Compton  electron  is 
generally  relativistic;  therefore,  the  energy  is  expressed  as 

m^c"’ 

E  =  KE  +  n,c'’  =  , - -  (2.6' 

\  1  -  V 

where  KE  -kin'  tic  energy 

'I 

-rest  mass  energy  of  an  electron 


KE  =  E 


Y  1+p 


where 


V 


-velocity  of  Compton  recoil  electron 


Solving  for  the  velocity 


V  = 


(2.7) 


The  primary  energy  loss  mechanism  for  Compton  electrons  is  coulomb 
collisions  with  molecular  electrons.  The  collisions  result  in  the 
ionization  of  air  molecules.  Experiments  show  that  one  electron-ion 
pair  is  produced  per  3^  electron  volts  of  energy  deposited  by  a  Compton 
electron.  A  IMeV  Compton  electron  will  create  30,000  electron-ion  pairs 
during  its  time  of  flight.  For  'nighly  relativistic  electrons,  these 
CO  1  lisions  ha  ve  1  it  t  le  ef  fee  t  on  the  C  ompton  e  lec  tron  velocity.  As  the 
Compton  electron  loses  energy,  the  deflections  from  coulomb  col  lisions 
become  large.  As  a  result,  the  "lompton  electron  has  a  forward  component 
of  velocity  close  to  the  speed  of  light  over  most  of  its  forward  range. 
The  Compton  electron  reaches  its  effective  range  when  the  deflections 
become  large  and  the  net  current  produced  goes  to  zero.  A  Compton 
electron  traveling  at  the  initial  velocity  over  the  effective  range  is  a 
convenient  approximation  when  modeling  the  Compton  currents. 

RETARDED  TIME 

Tne  gamma  ray  pulse  from  a  nuclear  detonation  can  be  visualized  as 
an  expanding  spherical  shell  with  a  thdekness  of  a  few  meters.  In  a 
standard  time  frame  of  reference,  the  arrival  time  of  the  pulse  wi  1  1  be 
a  function  of  distance  from  the  burst.  A  transformation  is  made  to  a 
frame  of  reference  where  a  given  phase  of  the  gamma  pulse  wi  1  1  arrive  at 
any  location  at  the  name  time.  The  modified  frame  of  reference  is 


called  retarded  time  (t).  Retarded  time  is  defined  as 

T  =  t  -  I  (2.8) 

c 

where  t  -standard  time 

r  -distance  from  burst 

c  -speed  of  light 

The  retarded  time  frame  is  convenient  because  the  phase  of  the  gamma 
pulse  does  not  shift  with  position. 

COMPTON  CURRENTS  AND  SECONDARY  ELECTRON  DENSITY 

Approximations  for  the  Compton  currents  and  secondary  electron 
density  were  developed  by  K  arzas  and  Latter  (Ref  9).  These  expressions 
are  based  on  two  assumptions.  First,  the  Compton  electrons  travel  at 
the  initial  velocity  over  the  entire  range.  Second,  the  Compton 
electron  production  rate  is  constant  over  the  range  of  a  Compton 
electron.  The  range  of  Compton  electrons  at  the  peak  of  the  deposition 
region  is  on  the  order  of  tens  of  meters.  Figures  2.1  thru  2.7  show 
tliat  Compton  electron  production  rates  vary  slowly  over  these  distances. 

The  expressions  are  best  understood  by  writing  them  down  as  they 
would  be  in  the  absence  of  a  magnetic  field. 

Compton  electron  density 
R/v^ 

n-(r,T)  =  g(r)  Y  /  f(T,T')  dx'  (2.9' 

^  I  o 

where  n|,(r,T)  -density  of  Compton  recoil  electrons 

Y^  -gamma  yield 

R  -Compton  electron  range 


-initial  Compton  electron  velocity 


f(T,T’)  -time  dependence  of  the  gamma  flux  accounting 
for  the  electron's  progression  through 
retarded  time 

f(T,T')  =  f[T  -  T'(l-Vj^/C)]  (2.10) 

The  gamma  ray  pulse  travels  at  the  speed  of  light.  A  phase  of  the  gamma 
pulse  will  have  a  fixed  value  of  retarded  time.  The  Compton  electrons 
are  traveling  at  Vq<C;  therefore,  the  electron  will  progress  through 
retarded  time  at  a  rate  of  (1-Vq/C).  During  the  Compton  electron's  time 
of  flight  (R/Vq),  the  electron  will  contribute  to  the  Compton  electron 
desnity  in  the  retarded  time  range  of  t  (birth)  to  t+ ( 1  - Vq/C)R/ v^. 
Equation  2.9  is  the  expression  for  the  Compton  electron  density  at  r  and 
T  where  the  integral  sums  the  contribution  from  earlier  retarded  time. 
Expressions  for  the  secondary  electron  production  rate  and  Compton 
current  are  based  on  the  Compton  electron  density. 

The  secondary  electron  production  rate  (ng)  from  one  Compton 
electron  is  assumed  to  be  constant  over  the  Compton  electrons  range. 
The  secondary  electron  production  rate  is  expressed  as 

V  R/v 

n„(r, t)  Q_  g(r)  /  f(T,T')  dx'  (2.11) 

^  R  o 

where  Q  -number  of  secondary  electrons  created  per  Compton 
electron 
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The  net  Compton  current  would  be  in  the  radial  direction. 

R/v 

Jp(r,T)  ^  -ev^g(r)e  /  °f{T,T')  di*  (2.12) 

O 

where  e  -average  cosine  of  the  initial  Compton  recoil 
electron  deflection  angle 
e  -charge  of  an  electron 

The  earth’s  geomagnetic  field  impacts  the  results  in  two  ways. 
First,  the  curvature  of  the  Compton  electron's  path  causes  an 
acceleration  of  its  movement  through  retarded  time.  The  expression  for 
the  time  dependence  of  the  gamma  flux  accounting  for  the  electron's 
progression  through  retarded  time  becomes 

f(T,T)  =  f[  r  -  (1  -  IcOs2(0))t'  +  ^  sin-(0)sin(,iiT' )]  (2.13) 

c  too 

where  0  -angle  between  the  earth's  geomagnetic  field  and  the 
direction  of  gamma  ray  motion 
w  -Larmor  frequency  for  the  Compton  electrons  in  the 
earth’s  geomagnetic  field 

Second,  transverse  currents  are  produced.  °olar  coordinates  with  the 
axis  along  the  earth's  geomagnetic  field  (Figure  2.S)  are  used  to 
describe  these  currents. 
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Figure  2.8  Polar  Coordinate  System 


Radial  Compton  Current 


R/v^ 

J„(r,T)  =  -eev„g(r)  /  °  f(T,T' )[C0s2(0) 

+  SIN“(0)COS(ut' )]dT'  (2.14) 

Theta  Compton  Current 

R/v^ 

jQ(r,T)  =  -ecv^g(r)  /  ^  f(  t  ,t' )STN(0)COS(O) 

•[COS(>.)T*  )  -  UdT'  (2.15) 

Phi  Compton  Current 

R/v^ 

J^(r,T)  =  -eev^g(r)  /  f(T,T' )SI\(0) 

•3IN(u)t' )dT’  (2.16) 

AIR  CHEMISTRY 

Air  chemistry  refers  to  the  processes  that  govern  the  secondary 
electron  and  ion  densities  (Ref  7).  The  processes  are  described  here  in 
order  of  dominance  in  early  time  for  high  altitude  EMP. 

IONIZATION,  Compton  electrons  ionized  oxygen  and  nitrogen  in  the 
air.  Ionization  is  the  source  for  secondary  electrons  and  positive 
ions.  On  the  average,  3^  electron  volts  of  energy  are  lost  by  the 
Compton  electron  per  electron-ion  pair  created. 

O2  +  ENERGY  —  0:^  +  e“ 

N2  +  ENERGY  —  N2  +  e" 

ELECTRON  ATTACHMENT.  Secondary  electrons  attach  to  oxygen 
molecules  to  form  a  negative  ion.  The  primary  mode  for  attachment 
requires  the  simultaneous  collision  of  the  secondary  electron  and  two 


oxygen  molecules: 


e  +  ^2  ^  ^2  *  ^2  ^  ^2 

Electron  attachment  removes  secondary  electrons  and  creates  negative 
ions.  The  rate  equations  describing  this  process  are 

d 

_  N  (  T  ,  r )  =  -K  ,  ( r ,  T  >  N„  ( r ,  T )  ( 2 . ; 

dx  c  I  e 

^N(r,T)  =  K,(r,T)?J{r,r)  [2.18) 

dx  -  I  e 

where  N^Cr,!)  -secondary  electron  density 

N  (r,T)  -negative  ion  density 

K^(r,T)  -electron  attachment  rate  constant 

The  electron  attachment  rate  constant  is  proportional  to  air  density 
squared  and  is  a  function  of  humidity  and  electric  field  strength.  For 
high  altitude  EMP,  the  humidity  is  close  to  zero. 

DISSOCIATIVE  RECOMBINATION.  Secondary  electrons  attach  to  positive 
ions.  The  electron  binding  energy  dissociates  the  molecule: 


e“  +  0^  —  0  +  0 
9~  +  N|  N  +  N 


Dissociative  recombination  removes  secondary  electrons  and  positive 
ions.  The  rate  equations  describing  this  process  are 


d 

d  1 


T  ) 


-k2Ng(r,T  )N_^(r,r  ) 


(2.19) 
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\ 

U 


—  N.(r,T)  =  -k5N„(r,T)N  (r,T)  (2.20) 

dT  ^  ^ 


where  K2  -dissociative  rate  constant 


The  dissociative  recombination  rate  is  approximately  equal  to 
2x10^  m^/sec  and  is  independent  of  air  density. 

NEUTRALIZATION.  Negative  ions  attach  to  positive  ions.  The 
process  requires  a  three  body  collision: 


o:t  +  Ot  +  M 


N-  +  N|  +  M 


^2  *  ^2 


N2  +  N2  +  M 


where  M  -third  molecule 

Neutralization  removes  positive  and  negative  ions.  The  rate  equations 
describing  the  process  are 


_  N.(r,T)  =  -ko(r)N  (r,T  )N  (r,T 
dx  J  * 


—  N  (r.x)  =  -l<o(r)N  (r,T)N  (r.i) 
dx  “  J  * 


(2.21  ) 


[2.22'' 


where  (r)-neutralization  rate  constant 

The  neutralization  rate  constant  is  proportional  to  air  density  and  is 
approximately  equal  to  1  ra^/sec  at  sea  level. 

The  approximate  values  for  the  dissociative  recombination  rate 
constant  and  the  neutralization  rate  constant  are  adequate  because  these 
processes  have  little  effect  on  early  time  secondary  electron  and  ion 
densities.  The  rate  equations  of  these  processes  are  combined  to  form  a 
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set  of  coupled  differential  equations  that  govern  the  secondary  electron 


and  ion  densities. 


d 

'dr 


Mg(r,T) 


ng(r,T)  -  (r,T)Ng(r,T) 
-k2Ng(r,T)N_^(r,T) 


(2.23) 


^  N_^(r,T)  =  ngCr.x)  -  k2Ng(r,T)M^(r,T) 

-k3(r)N^(r,T)N_(r,T)  (2.2^4) 

(r.x)  =  -k. (r,T)N-(r,T) 
dT  "  «  c 

-k2(r)N^(r,T)N_(r,T)  (2.25) 

AIR  CONDUCTIVITY 

The  secondary  electrons  and  ions  drift  in  response  to  the  electric 
field  (E).  Air  conductivity  (o)  relates  the  displacement  current 
density  (J)  to  the  electric  field  strength: 

3  s  0  ^  (2.26) 

Mobility  (y)  is  used  to  relate  the  drift  velocities  to  the  electric 
field  strength.  The  mobility  for  ions  is  approximately  equal  to 
2.5x10“^  m^/(volt  sec)  for  STP  (standard  temperature  and  pressure)  air 
and  is  inversely  proportional  to  air  density.  The  electron  mobility  is 
larger  than  the  ion  mobility  and  is  a  function  of  field  strength  and 
humidity  as  well  as  being  inversely  proportional  to  air  density.  The 
conductivity  expressed  in  terms  of  charge  density  and  mobility  is 

o(r,T)  r  [Ng(r,T)  yg(r,T)  +  {N^(r,T)  +  N^(r ,  t)  )u  j^(  r )  ]q  (2.27) 
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where  Mg(r,T),Uj^(r)  -electron  and  ion  mobilities 

q  -charge  of  an  electron 

Air  conductivity  is  primarily  dependent  on  electron  mobility  and 
density . 


ELECTROMAGNETIC  FIELDS 

The  electromagnetic  fields  are  described  by  Maxwell's  equations. 
Maxwell's  equations  expressed  in  rationalized  MKS  units  are 


(2.28) 


V  X  E 


i5 

at 


(2.29) 


V  •  §  =  0 


(2.30) 


7x5 


if 

at 


where 


E  -electric  field  vector 

B  -magnetic  field  vector 

p  -charge  density 

-permittivity  of  free  space 
y  -permeability  of  free  space 

o 

j  -current  density  vector 

c  -speed  of  light 


(2.31) 
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Maxwell's  equations  are  expressed  in  polar  coordinates  and  retarded  time 


by  making  the  following  transformations: 

3 _ ^  3_ 

3t  3t 


where 


e 

r 


-unit  vector  in  the  radial  direction 


Maxwell's  equations  expressed  in  polar  coordinates  and  retarded  time  are 
written  as 


V  • 


1  3  ♦ 
c  3t^ 


0 

£ 


(2.32) 


V  X  E  -  e 


^  3  ?  9  ^ 

*  c  3t^  '  3t® 


(2.33) 


V  •  S  -  e 

r 


1  3  + 

-  r  0 

c  3t 


(2.3^*) 


V  X  B 


e 

r 


1  3  - 

X  -  ^B 

c  3t 


j 


(2.35) 


Karzas  and  Latter  (Ref  9)  derived  differential  equations  that  describe 
the  transverse  and  radial  electric  fields.  Terras  involving  spatial 
derivatives  are  neglected  when  compared  to  derivatives  with  respect  to 
T.  This  is  known  as  the  high  frequency  approximation.  The  differential 
equations  resulting  from  the  high  frequency  approximation  are  written  as 
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I 


(2.36) 


(2.37) 


(2.38) 


where 

Eq,  Ej|j,  -components  of  electric  field  vector 

Jq,  JjJj,  Jj,  -components  of  current  density  vector 

The  current  density  is  the  vector  sum  of  the  Compton  current  density  and 
the  conduction  current  density.  The  conduction  current  is  the  current 
from  electrons  and  ions  drifting  in  the  local  electric  field. 

The  differential  equation  for  the  radial  electric  field  (EQ  2.36) 
describes  the  time  variation  of  the  radial  electric  field  in  terms  of 
the  local  current  density.  The  differential  equations  for  the 
transverse  electric  field  (EQS  2.37  and  2.38)  describe  the  spatial 
variation  of  the  transverse  electric  fields  in  terms  of  the  local 
current  density  and  conductivity.  The  transverse  electric  fields  are 
generated  in  the  absorption  region  and  are  able  to  propagate  to  the 
ground. 
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III.  NUMERICAL  METHODS 


Chapter  III  discusses  the  numerical  methods  used  in  EMPFRE  to  model 
the  processes  discussed  in  Chapter  II.  The  approach  used  to  solve  for 
the  electric  fields  was  developed  by  J.  H.  Erkkila  (Ref  5).  Many  of  the 
approximations  made  to  model  the  phenomenon  were  proposed  by  W.  J. 
Karzas  and  R.  Latter  (Ref  9).  The  details  of  the  methods  used  are 
presented  in  the  following  sections. 

BURST  GEOMETRY 

The  geometry  and  coordinate  system  used  for  EMPFRE  are  shown  in 
Figure  3.1.  ZRMIN  and  ZRMAX  are  the  upper  and  lower  bounds  of  the  gamma 
ray  absorption  region.  RMIN  and  RMAX  are  the  distances  from  the  burst 
point  to  ZRMIN  ani  ZRMAX  along  R.  EMPFRE  calculates  the  electric  fields 
along  R  from  RMIN  to  RMAX.  The  transverse  componenets  of  the  electric 
field  propagate  from  RMAX  to  the  target  with  a  loss  due  to  spherical 
divergence , 

If  the  height  of  burst  is  less  than  the  default  value  for  ZRMIN  (50 
kilometers),  ZRMIN  is  set  equal  to  the  height  of  burst  and  ZRMAX  is  set 
equal  to  the  target  altitude.  If  the  target  altitude  is  greater  than 
the  default  value  for  ZRMAX  (10  kilometers),  ZRMAX  is  set  equal  to  the 
target  altitude.  If  ZRMAX  is  set  to  the  target  altitude,  the  radial 
component  of  the  electric  field  is  included. 

Polar  coordinates  centered  at  the  burst  point  are  used  to  compute 
the  electric  fields.  The  polar  axis  is  parallel  to  the  earth's  magnetic 
field. 
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Pol  ar 
rty  1  s 


The  target  location  is  designated  by  the  rectangular  coordinates  x, 
y  and  z  in  meters.  The  earth  is  assumed  to  be  flat  and  z  is  the  target 
altitude.  The  x-axis  points  to  magnetic  west.  The  y-axis  points  to 
magnetic  south. 

BANGLE  (the  geomagnetic  field  dip  angle)  is  specified  in  degrees. 

The  angle  theta  (9)  is  found  by  taking  the  dot  product  of  a  unit 

/ 

vector  in  the  direction  of  the  earth's  magnetic  field  (8)  and  a  unit 
vector  in  the  radial  direction  (r): 

GrARC  C0S(B*r)  (3.I) 

where  8=C0S(BANGLE)  j  +  SIN(BANGLE)  J< 
r=x  1  +  y  j  +(H0B  -z)  ft 

GAMMA  RAY  YIELD  FUNCTION 

EMPFrtE  is  designed  to  analyze  the  electro  .'nagnetic  pulse  from  an 
arbitrary  mono-energetic  gamma  pulse.  The  total  gamma  ray  yield 
(kilotons),  gamma  ray  energy  and  a  function  (F(t))  that  is  proportional 
to  the  gamma  ray  yield  as  a  function  of  time  are  required  as  input 
parameters.  The  yield  in  kilotons  is  converted  to  MeV.  The  gamma  r'r.-j 
energy  is  used  to  calculate  the  yield  as  number  of  gamma  rays  emitted. 
Finally,  a  proportion,!  1  i  ty  constant  is  foun.1  to  convert  the  value  of  tne 
function  to  the  gamma  ray  emission  rate.  The  proportionality  constant 
is  derived  by  integrating  the  function  using  the  trapezoidal  rule: 

O 

where  F(r)  -gamma  function 

ti  -retarded  time  step 


pa 


T=r:  •  At 


The  integral  of  the  chosen  function  multiplied  by  the  proportionality 
constant  must  be  equal  to  the  total  gamma  ray  yield; 

/“c*F(T)dT=YIELD  (3.3) 

O 

The  constant  can  be  taken  out  of  the  integral  and  solve::  for  witn  the 
integral  approximated  by  the  trapezoidal  rule. 


YIELD 

At-[F(t, )/2+F(t2) ...♦F(t^)/2] 


(3.^) 


The  function  F(  t)  can  be  in  the  form  of  a  mathematical  expression  or  a 


set  of  data  points  where  interpolation  is  used. 


INITIAL  ELECTRON  VELOCITY 

The  intitial  relativistic  energy  of  the  Comptor.  electron  is  given 
by 

ErKE^+m^c^  (3.5' 

Erra^c^/C (3.6 j 
where  -initial  kinetic  energy 

Vq  -initial  velocity 

Solving  for  the  initial  velocity,  the  equation  is 

VQ=2.99E8-SQRT(1-(0.5n /(KEq+0.511 ))^)  (3.7' 

where  -initial  kinetic  energy  (MeV) 

-initial  velocity  (meters/second) 

MEAN  FREE  PATH 

The  absorption  coefficients  (w )  for  Y-ray  interactions  in  air  are 

shown  in  Figure  3.2.  The  mean  free  path  is  given  by 

MFP  =  _1  (3.=^) 

UP 


where 


p  -the  density  of  air 


II  Air 


The  mean  free  path  for  Compton  interaction  and  total  absorption  is  found 
using  an  empirical  approximation  to  Figure  3.2.  The  equations  used  to 
approximate  the  mean  free  path  is  STP  air  as  a  function  of  energy  are 


For  0.2  MeV  <  S  <  0.5  MeV 

MFP^  =  MFPj_=l  /[  0.0035+ (In  (E)+1 .609) /3271  ]  (3.9  ) 

For  0.5  MeV  <  E  <  1.5  MeV 

MFP^=MFPj.=  1/[  0.0038-(ln(E)+0.5931  )/1927]  (3.10} 

For  1.5  MeV  <  E  <  10.  MeV 

MFP^=1/[0.0033+(0.9055-ln(E 'O/gSS]  (3-11  ) 

t-IFP^=l  /  [ 0.0033+ (0.^055-ln(E) ;  ]  (3.12) 


where  E  -ga:nna  ray  energy  (Me"' 

I'^FPj^  -nean  free  path  for  Compton  interaction 
MFPj.  -mean  free  path  for  all  interactions 

The  MF?  is  inversely  proportional  to  air  density,  therefore;  MF?  as  a 
function  of  altitude  is 

MFP(ALT)  =  MFP3-:,  EXP(ALTir.JD£/7000  Meters^  (3.13'' 

The  MFP  for  Compton  interaction  and  total  absorption  (where  different) 
as  a  function  of  ganma-ray  energy  in  STP  air  is  generated  from  an 
empirical  fit  of  Figure  3*2.  Values  for  the  MF'’  in  STP  air  are  provided 


in  Table  3.1 


Gamma 

Energy 

(MeV) 

Compton 

MFP 

(m) 

Total 

MFP 

(m) 

Compton 

Electron 

Range 

(m) 

Initial 

Electron 

Velocity 

(v/c) 

0.2 

285.7 

0.021 

0.37 

i 

0.3 

275.9 

0.073 

0.50 

O.U 

269.4 

0.16 

0.59 

0.5 

264.6 

0.28 

0.67 

0.75 

282.5 

0.68 

0.78 

1.0 

295.0 

1 .21 

0.85 

2.0 

329.2 

319.7 

4.06 

0.95 

3.0 

380.9 

351.8 

7.49 

0.98 

4.0 

428.8 

378.7 

11.2 

0.99 

5.0 

475.0 

402.6 

15.0 

0.992 

590.7 

454.8 

24.9 

0.997 

714.2 

500.8 

34.0 

0.998 

Table  3.1  Data  Used  in  EMPFRE 


COMPTON  DYNAMICS 


Values  for  the  median  Compton  electron  energy,  Compton  electron 
deflection  angle  and  gamma  ray  deflection  angle  were  found  by  using  the 
K line-Nishina  cross  sections  for  Compton  scattering.  First,  the  median 
Compton  electron  deflection  angle  is  determined.  The  corresponding 
Compton  electron  energy  and  gamma  ray  deflection  angle  are  calculated. 
COMP  is  a  program  written  by  LTC  J.  Erkkila  that  performs  this  task. 
Empirical  fits  of  the  results  are  used  in  EMPFSE  and  are  given  below: 


Compton  Electron  Energy 
For  E.  <  4  MeV 

E=Ev  .  (0. l95+(ln(E. )+1 .609)/5.992) 
For  MeV 

S=E  •  (0.b95+(ln(E|,)-1.386).'9.163) 

Compton  Electron  Deflection  Angle 
For  0.2  MeV  <  E^  <  0.5  MeV 

0=1.2488-0. 18« (ln(EY)+1 .6094) 

For  0.5  MeV  <  <  1.5  MeV 

O:' ,0383-0. 1422- (ln(E, )+0. 6931) 

For  1.5  MeV  <  E,^  <  10.  MeV 

0=0.91  63  +  0.137  K0.4055-ln(E|) 


Gamma  Ray  Deflection  Angle 
For  0.2  MeV  <  E^  <  0.5  Mev 


;|,  =  U.YO04-0.U:^‘jc/ 


^  All  n, 


Y 


;+ 1  .  oJ:# ! 


(3.14) 


(3.15) 


(3.16) 


(3.17) 


(3.18) 


.  ly  / 


I 


For  0.5  MeV  <  E^<  10.  MeV 

<j.=0. 698 1-0. 6861  •  (ln(E^)+0.6931  )  (3.20) 

where  E  -Compton  electron  energy  (MeV) 

-gamma  energy  (MeV) 

0  -Compton  electron  deflection  (radians) 

?  -gamma  deflection  angle  (radians) 

Values  of  Compton  electron  energy,  Compton  electron  deflection  angle  and 
gamma  ray  deflection  angle  from  the  empirical  fits  are  provided  in  Table 
3.2  along  with  the  results  from  COM®. 


COS  0 
[Ref] 
(Rad) 

Error 

(i) 

COS  4. 
[Ref] 
(Rad) 

Error 

(It) 

Energy 

[Ref] 

(MeV) 

0.3165 

[0.3160] 

0.16 

0.7071 

[0.7082] 

0.15 

0.0390 

[0.0422] 

0. 38^18 
[0.3849] 

0.03 

0.7338 

[0.7267] 

0.98 

0.0788 

[0.0796] 

0.4320 

[O.4313] 

0.16 

0.7522 

[0.7671] 

0.68 

0.1242 

[0.1232] 

0.4679 

[0.4650] 

0.62 

0.7660 

[0.7671] 

0. 14 

0.1739 
[0. 1718] 

0.5181 

[O.5213] 

0.61 

0.8179 

[0.8107] 

0.89 

0.3116 

[0.3095] 

0.5527 

[0.5576] 

0.88 

0.8511 

[0.8443] 

0.81 

0.4635 

[0.4640] 

0.6396 

[0.6370] 

0.41 

0.9182 

[0.9180] 

0.02 

1.158 

[1.174] 

0.6813 

[0.6797] 

0.24 

0.9485 

[0.9487] 

0.02 

1.941 

[1.959] 

0.7096 

[0.7086] 

0.14 

0.9659 

[0.9644] 

0.16 

2.780 

[2.781] 

0.7308 

[0.7300] 

0.11 

0.9769 

[0.9735] 

0.35 

3.597 

[3.677] 

0.7676 

[0.7668] 

0.10 

0.9915 

[0.9849] 

0.67 

5.727 

[5.804] 

Error 

(J) 

7.5 

1.0 

0.81 

1.2 

0.67 

0.11 

1.U 

0.92 

0.00^4 

2.2 

1.3 


0.7923 

[0.7909] 


0.18 


0.9976 

[0.9900]  0.77 


7.950 

[8.036] 


1.1 


ELECTRON  RANGE 


Range  equations  are  from  Radiological  Healtn 

For  0.01  MeV  <  E  <  2.5  MeV 

R=14 12. E'' (1  .265-0. 095^1 -InCE)) 

For  E  >2.5  MeV 
R=530.E-106 

where  R  -range  (mg/cm) 

E  -energy  (MeV) 

The  range  in  STP  air  is  found  by  dividing  the 
density  of  ST?  air  (P2-pp:1.293  mg/cm^).  Conv 
correcting  for  altitude,  the  final  equations  are 

For  0.01  MeV  <  E  <  2.5  MeV 

R=3. 19*E'(1.265-0.095il-ln(E))  .EXP(ALT/7000) 

For  E  >  2.5  MeV 

R=4 . 10. E-0. 82. EXP (ALT/7000) 

where  R  -range  (meters) 

ALT  -altitude  (meters) 

Electron  range  as  a  function  of  energy  in  ST?  ai 


Handbook  'Ref  16 ) . 


(3.22) 


range  (above)  by  the 
ertir.g  to  meters  and 

(3.23) 

(3.2il) 


r  is  provided  in  Table 


3.2.  The  electron  range  includes  the  effects  of  straggling;  therefore 
the  range  is  the  forward  component  of  the  electron's  path: 


N 

R=  I  PATH  C3.25) 

n=  1 

where  PATH  -distance  traveled  between  collisions  n-1  and  n 

N  -total  number  of  collisions  to  themalize  the  electron 
Q  -deflection  angle  from  original  path 

The  average  cos(Q^)  term  approaches  1  for  high  energy  electrons  early  in 
TOP  (time  of  flight)  and  goes  to  zero  for  low  energy  electrons  after 
many  collisions. 

COMPTON  ELECTRON  TIME  ^  FLIGHT 

The  Compton  electron  time  of  flight  Is  approximated  by  assuming  the 
Compton  electron  travels  at  its  initial  velocity  over  its  entire  range. 
Therefore,  the  Compton  electron  time  of  flight  is  written  as 

TOF=RANGE(R)/v^  (3.26' 

o 

TOFrRANGE^  • EXP ( ALT/7000 ) /v^  (3.27) 

where  HANGE(R)  -Compton  electron  range  at  fi 
-initial  velocity 
RANGEq  -range  in  oTP  air  (meters) 

ALT  -altitude  (meters) 

Assuming  a  constant  velocity  over  the  entire  range  is  justified  for 
relativistic  Compton  electrons  at  high  altitudes  because  energy  loss  in 
early  TOF  primarily  reduces  the  electron’s  mass.  Later  in  the  electron 
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lifetime,  the  electron  will  have  lagged  behind  the  pulse  in  retarded 
time  primarily  due  to  its  curved  trajectory.  Electrons  lagging  behind 
the  gamma  ray  pulse  become  unimportant  due  to  the  increase  in 
conductivity. 

COMPTON  ELECTRON  PRODUCTION  RATE 

The  Compton  electron  production  rate  is  calculated  with  the 
assumption  that  scattered  gamma  rays  are  ignored.  The  gamna  ray  flux  is 
expressed  as 

R 

EXP[-/drVA^(r'  )] 

F  (R,t)  =  S^(t)  - 2 _  (3.2B' 

^  ^  Mr,R2 

where  X(r)  -gamma  ray  MFP  at  r  (met-ics' 

R  -position  along  slant  range  'meters; 
s^(t)  -gamma  ray  emission  rate  at  t  isec"^) 

Assuming  an  exponential  atmosphere,  the  MFP  as  a  function  of  r  can  be 
related  to  the  MFP  in  STP  air  as 


MFP(R)=MFPg-EXP 


■ H03-R-C0S(A)‘ 

— 7m — 


(3.29) 


where 


MFPq  -mean  free  path  for  STP  air  (meters) 

HOB  -height  of  burst  (meters) 

A  -angle  between  vertical  and  slant  range 


i 


I 


The  Compton  electron  production  rate  is  found  by 


N^,(R,t)  =  S^(t) 


EXP[-dr/MFP^.(R)] 


^♦TtR^-MFPgCR) 


(m'^sec'^ ) 


(3. 30; 


where  MFPj.(r)  -mean  free  path  for  total  absorption 

MPP^Cr)  -mean  free  path  for  Compton  interaction 

For  an  exponential  atmosphere,  the  formula  for  Compton  electron 
production  rate  as  a  function  of  R  is  expressed  as 


EXP 


N^(R,t)=S^(t) 


■  -7000 

/ 

■  R-C0S(A)-H0B  ■ 

■  -HOB' 

\ 

(EXP 

-EXP 

1 

COS(A)  •  MFP^ 

7000 

_  7000  _ 

j 

4tir2.mFP.  -  EXP 


HOB-R-COS(A)' 


7000 


- '  (3.31  ) 


where  ^^c  free  path  for  STP  air  'meters) 

The  gamma  rays  are  assumed  to  be  in  a  collimated  beam.  A  scattered 
gamma  ray  changes  direction  and  lags  behind  the  pulse.  Compton 
electrons  produced  by  these  gamma  rays  will  not  contribute  to  the  EMP 
because  of  the  increase  in  conductivity  at  late  times.  Therefore, 
scattering  removes  the  gamma  ray  from  the  pulse. 
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COMPTON  CURRENT  AND  SECONDARY  ELECTRON  SOURCE 


Integral  equations  for  the  Compton  current  density  and  the 
secondary  electron  source  density  developed  by  Karzas  and  Latter  are 
solved  numerically  (Ref  9). 

Secondary  Electron  Production  Rate 


TOP 


NgCr.T)^  Q  G(r)  /  f(T,T')dT' 


(3.32) 


Radial  Compton  Current  Density 


TOE 


-evQG(r)  /  f(T,T')[COs2(0)+SIN2(G)COS{wT’)]d' 


(3.33) 


Theta  Compton  Current  Density 


TOE 


Jg(r,T)a-evQG(r)  /  f  (t  ,t' )  [SIN(0)C0S(0)  (COS(a)T' )-1 )  Idi’ 


(3.3^4) 


Phi  Compton  Current  Density 


TOP 


J^(r,T)%-eVgG(r)  f(T,  T)tSIN(e)SIN(wf  )3dT' 


(3.35) 


where  TOP  -Compton  electron  time  of  flight 


f(T,T')  -time  variation  of  gamma  emission  (sec“^) 


o  o  SIN(  uf) 

f(T,T')  =  f[T-(1-P*COS^(0))*T'+3*SIN^(O)  - ] 


(3.36) 


T  -retarded  tiffle(sec) 
r  -position  along  the  slant  range  (ro) 

Q  -number  of  secondary  electrons  produced  per 


Compton  electron 

Vq  -initial  Compton  electron  velocity  (m/sec) 

R(r)  -Compton  electron  range  (m) 

G(r)  -Compton  electron  source  density  per 
yield  (m“^) 

e  -electron  charge  (coulombs) 

0)  -Larmor  frequency  for  a  Compton  electron  (sec“^) 

0  -the  angle  between  the  slant  range  and  the  geomagnetic  field 
6  -initial  velocity  (fraction  of  c) 

The  integral  equations  were  solved  numerically.  The  equations  used 
to  solve  for  Compton  current  density  and  the  secondary  electron  source 
density  are  listed  below: 

Secondary  Electron  Production  Rate 

v  n 

NQ(r,T)a:Q  G(r)[  f(T,Tj^')  Ax'  (3.37) 

Radial  Compton  Current  Density 

n 

Jr(r,T)a-ev'G(r)  I  f  (t  ,1^,' )  [COs2(0)+SIN2(0)COS(a)Tj,’ )  ]At  ’  (3.38) 

n=  1 

Theta  Compton  Current  Density 

n 

J^^r,T)s: -ev'G(R)  I  f  (  )  [SIN(0)COS(0)  (C0S(u)t^»  )-1 )  ]At'  (3.39) 

n=  1 


I 
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Phi  Compton  Current  Density 

n 

J,(r,T)s:-ev’a(r)  V  f(T  )SIN(0)SIN(u)t^' )At'  (3.40) 

n=1 

where  At'  -retarded  time  step  size  (0.05  Shakes) 

T^'  =n*AT' 

N  =R(r)/v^/AT' 

v'  -component  of  the  initial  velocity 
in  the  radial  direction 

These  computations  are  performed  in  subroutine  current.  The 

retarded  time  step  used  to  approximate  the  integrals  is  0.05  shakes. 

The  gamma  output  of  the  device  will  change  slowly  over  this  time  for 
most  gamma  histories.  This  subroutine  drives  the  computation  time. 
The  program  uses  about  seconds  of  cpu  time  per  electric  field 

evaluation  on  the  CDC  CYBER  750  at  this  time  increment.  Decreasing  the 

step-size  becomes  very  expensive. 
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AIR  CHEMISTRY 


The  set  of  differential  equations  which  govern  the  electron  and  ion 
densities  as  a  function  of  time  are  written  as 
Electron  Density 


d 

— N 

dT 


e(r,T) 


S(r,T)-K^ (r)»Ng(r,T)-K2 


Ne(r,T)*N^(r,T) 


G.^1) 


Positive  Ion  Density 


— N^(r,T)  r  S(r,T)-K2*Ng(r,T)»N^(r,T)-K2(r)*N^(r,T)»N_(r,T) 


(3.^2) 


Negative  Ion  Density 


^_(r,T; 

dT 


K-i  (r)  •Ng(r,T)-K2(r)  •N^(r,T)-N_(r,T) 


(3.^3) 


where 


Ng(r,T) 

N^(r,T) 

N_(r,T) 

S(r,T) 

(r) 

^2 

Kjir) 


-electron  density  (m"^) 

-positive  ion  density  (m'^) 

-negative  ion  density  (m~') 

-secondary  electron  source  density  (m“3*sec“^) 
-attachment  rate  stant  (sec“^) 

(proportional  to  air  density  squared) 
-dissociative  recombination  rate  constant  (m^/sec) 
-neutralization  rate  constant  (m^/sec) 
(proportional  to  air  density) 


The  secondary  electron  source  density  is  calculated  at  intervals 
equal  to  the  electric  field  retarded  time  step.  Its  value  is  assumed  to 
be  a  linear  function  of  time  between  evaluations.  The  first  derivatives 


with  respect  to  retarded  time  are  replaced  by  a  forward  difference 


approximation.  The  numerical  method  used  to  solve  the  air  chemistry 
differential  equations  are  written  as 
Electron  Density 

Ne(r,Tn^l)  =  Ng( r , T^)+[S(r , (r) • Ng(r, 

-K2*Ng(r,T^)*N_^(r,T^)]*AT 

Positive  Ion  Density 

N+(r’.Tn+i^  =  N^(r',T^)+[S(r,T^)-K2*Ng(r,T^)*N^(r,Tj^) 

-K2(r)"N_j^(r,Tj^)*N_(r,T|^)3*AT  (3*^5) 

Negative  Ion  Density 

=  N_(r,T^)  +  [Ki  (r)«Ng(r,Tjj) 

-K2(r)‘N^(r,Tj^)*N_(r,Tj^)  ]*At 

whe-e  At  -retarded  time  increment 
"n 

The  retarded  time  increment  is  one  hundredth  the  size  of  the  increment 
used  for  the  electric  field. 

The  rate  constant  for  electron  attachment  is  a  function  of  the  air 
density,  the  electric  field  strength  and  the  water  vapor.  The  empirical 
fit  for  the  rate  constant  was  developed  by  Longley  and  Longmire  (Ref  11) 
and  employs  cgs  units.  The  electric  field  ( volts/meterl  is  converted  to 
electrostatic  units  per  atmosphere.  The  equations  were  modified  to  take 
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advantage  of  an  exponential  atmosphere.  The  equations  for  the  empirical 


fit  are  written  as 


e  =  3.331x10"5*E*EXP(ALT/7000) 


where  E  -electric  field  (volts/meter) 

ALT  -altitude  (meters) 

£  -electric  field  (esu/atmosphere) 


For  e  <  e . 


Eg  =  /(1+2.457- 


For  <  e  <  £2 


£q  =  [(0.1l85*p^-^^®+e)°-5-0.3^'42*p°-®3^]2 


For  £2  <e 


EqS  -1.195*p°*®3^ 


(3.'^7) 


(3.^8) 


(3.'*9) 


(3.50) 


where 


£■]  =  0.07853*  ( 1+2. 457*p®*^^^  )  (esu/atraosphere) 
£2  =  3.015+1 • 195*p®’®^^  (esu/atmosphere) 


p  -percent  of  water  vapor 


Rate  constant  for  electron  attachment  (sec“^) 


(3.51) 

(3.52) 


K  =  (1-0.01  p)[1+0.3'1^-p)  A3+A2] 


(3.53) 


where  A2  -two  body  attachment  rate  (sec~  ) 


A2  - 1. 22x10®*EXP( -ALT/7000) •EXP( -21. 15/eo) 


(3.5M: 


-three  body  attachment  rate  (sec“^) 


1x10®-EXP(-2  ALT/7000)- (0.62+800  e|) 

1  +  1x103.e2.[eQ*(U0.03*e2)]0-33 


(3.55) 


The  rate  constant  for  dissociative  recombination  is  independent  of 
air  density  and  has  a  value 


K2  =  2x10“^^  (m^/sec) 


(3.56) 


The  rate  constant  for  neutra lization  is  proportional  to  air 
density.  The  equation  for  the  rate  constant  is  written  as 


=  EXP (-ALT/7000)  (m^/sec) 


(3.57) 


The  time  scale  that  EMPFRE  works  with  is  of  the  order  of  1x10“'^ 
seconds.  For  this  time  scale,  electron  attachment  is  the  key  process; 
therefore,  more  effort  is  concerned  with  modeling  this  process. 
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ELECTRON  MOBILITY  AND  CONDUCTIVITY 


The'electron  mobility  is  a  function  of  air  density,  electric  field 
strength  and  water  vapor.  The  empirical  fit  for  electron  mobility  was 
developed  by  Longley  and  Longmire  (Ref  16).  The  electron  mobility  is 
written  as 

3.331x10-'7-  P3 

Mg  =  -  (3-58) 

1-0.01  P+0.01-P-R 

where 

16.8+Cg 

Mg  =  1x1&°‘EXP(ALT/7000)*  _  (3.59) 

0. 63+26. 7*eg 

R  =  1. 55+210/(1+11. S-Cq  +7. 2*eQ2)  (3-60) 

p  -percent  of  water  vapor 

The  parameter  (e^)  is  defined  in  the  section  on  air  chemistry.  Table 
3.3  provides  values  of  the  electron  mobility  as  a  function  of  altitude 
and  electric  field  strength  with  no  water  vapor. 

The  positive  and  negative  ion  mobility  is  approximated  by 

Mi  =  2.5x10~^*EXP(ALT/7000)  (m^/volt/sec)  (3.61) 

The  conductivity  (amps/m/volt)  is  written  as 

6(r,T)=[Ng(r,T)*Mg(r,T)+(N^(r,T)+N^(r,T)).Mi(r)] 'q  (3.62) 

where  Ng(r,T)  -secondary  electron  density  (m“^) 

N^(r,T)  -positive  ion  density  (m“3) 

N_(r,T)  -negative  ion  density  (m"^) 

q  -electron  charge  (1.6x10“^^  coulombs) 


For  early  time,  the  conductivity  is  primarily  a  function  of 


secondary  electron  density  and  mobility;  therefore,  more  effort  is  used 
to  model  electron  mobility  than  for  ion  mobility. 

SOLVING  FOR  THE  ELECTRIC  FIELDS 

Three  differential  equations  are  solved  to  find  the  components  of 
the  electric  field: 

For  the  radial  electric  field  (E^,) 

6Ej.  Jp+a*Ej. 

—  =  -  -  (3.63) 

^o 


For  the  theta  electric  field  (E^) 


fir 


Wo'C 

-  -(jg+U'Eg) 

2 


(3.64) 


For  the  phi  electric  field  (E^) 


Po*c 


) 


(3.65) 


where  t  -retarded  time  in  seconds 

J  -radial  Compton  current  density 
Eq  -permittivity  of  vacuum 
(8.854x10"''2  farad/m) 
o  -conductivity 

Pq  -permeability  of  vacuum  (HirxIO"^  Henry/m) 
c  -speed  of  light  (2.99x10®m/sec) 


t 


— 

E-field 

ALTITUDE 

! 

(V/m) 

o 

B 

1  Kra 

10  Km 

20  Km 

40  Km 

60  Km  j 

0 

2.389 

2.756 

9.968 

41.59 

724.2 

12609 

1 

0.5 

2.388 

2.754 

9.950 

41.29 

644.8 

4981.  ! 

I 

5 

2.379 

2.742 

9.796 

38.80 

365.2 

1460.  i 

1 

50 

2.293 

2.629 

8.539 

25.74 

114.1 

464.4  1 

1 

500  , 

1.735 

1.929 

4.385 

8.899 

33.80 

t 

272.1  1 

i 

1 

1 

5K 

0.687 

0.735 

1.315 

2.539 

16.71 

247.8 

50K 

0.195 

0.208 

0.404 

1.062 

14.35 

145.3  j 

1 

Table  3.3  Electron  Mobility 
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These  differential  equations  were  solved  numerically  with  a  fourth 
order  Runge-Kutta  technique.  The  differential  equations  are  of  the  form 

6E 

—  =  f(x,E)  (3.66) 

6x 

The  Runge-Kutta  method  solves  equation  3.66  for  E^  at  Xq+Ax  when 
the  value  of  E^  is  known  at  x^  (Ref  5).  The  fourth  order  Runge-Kutta 
technique  generates  4  estimates  of  the  change  in  E  over  the  Interval  Ax: 


K1=Ax*f(xQ.Eo)  (3.67) 

K2=Ax*f(xQ+Ax/2,EQ+K1/2)  (3.68) 

K3=:Ax*f(xQ+*x/2,E^+K2/2)  (3.69) 

K4  =  Ax*f(xQ-i-Ax,E£)+K3)  (3*70) 


The  value  of  E^=E(Xq+  Ax'  is  approximated  by 

E^  =  Eq+(KU2  •(K2+K3)+K4)/6  (3.71) 

The  Runge-Kutta  method  is  convenient  for  two  reasons.  First,  the 
method  requires  one  value  of  the  electric  field  to  start  the 
computation.  Second,  the  independent  variable  increment  can  be  varied 
during  the  computation.  The  drawback  is  that  the  error  is  difficult  to 
estimate  (Ref  13:74).  When  the  second  derivative  of  the  electric  field 
with  respect  to  the  independent  variable  (x)  becomes  large,  the 
independent  variable  Increment  must  be  reduced  to  avoid  instabilities 


and  reduce  error 


Two  scenarios  lead  to  instabilities  in  the  Runge-Kutta  technique 


for  this  application. 

1)  Solving  for  the  radial  electric  field  differential 
equation  in  early  time,  the  second  derivative  with  respect  to  time 
will  be  positive  and  can  become  large.  An  instability  occurs  when 
the  initial  time  increment  causes  a  large  error  in  K1  (too  small). 

The  value  of  K2  (too  large)  over  corrects  for  K1.  The  value  of  K3 
(too  small)  over  corrects  for  K2  and  becomes  even  smaller  than  K1. 
Finally  (too  large)  over  corrects  for  K3  and  becomes  larger 
than  K2.  The  estimates  of  the  increase  in  the  radial  electric 
field  (K1-K4)  will  oscillate  when  the  second  derivative  is  not 
equal  to  zero.  These  oscillations  are  reduced  by  decreasing  the 
increment  when  the  second  derivative  becomes  large. 

2)  Solving  for  the  transverse  electric  field  as  saturation 
is  reached,  the  second  derivative  with  respect  to  position  will  be 
negative  and  can  become  large.  An  instability  occurs  when  the 
initial  position  increment  causes  a  large  error  in  K1  (too  large). 

The  value  of  K2  (too  small)  over  corrects  for  Kl.  The  value  of  K3 
(too  large)  over  corrects  for  K2  and  becomes  even  larger  than  K1. 
Finally  K4  (too  small)  over  corrects  for  K3  and  becomes  smaller 
than  K2.  The  estimates  of  the  increase  in  the  transverse  electric 
field  (K1-K4)  will  oscillate  when  the  second  derivative  is  not 
equal  to  zero.  These  oscillations  are  reduced  by  decreasing  the 
position  increment  when  the  second  derivative  becomes  large. 

Initial  development  of  EMPFRE  used  as  a  fixed  independent  variable 
increment  to  solve  the  differential  equations.  This  method  was  unstable 
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for  gamma  yields  that  drove  the  electric  field  to  saturation  because  the 
second  derivative  with  respect  to  the  independent  variable  becomes  large 
at  that  point.  Instabilities  are  forecasted  by  monitoring  the  behavior 
of  K1  thru  K4. 

When  the  oscillations  grow,  the  independent  variable  increment  is 
reduced  and  the  procedure  is  repeated.  Three  tests  were  used  to 
forecast  instabilities.  When  one  of  the  following  conditions  were  met, 
the  independent  variable  increment  was  halved: 


1) 

Ki2  > 

2 

K2^ 

(3.72) 

2) 

K2^  > 

2 

Ki2 

(3.73) 

3) 

K1  K4 

< 

0 

(3.74) 

An  instability  observed  without  these  checks  is  modeled  with  a 
program  designed  to  study  the  behavior  of  K1-K4  and  the  results  are 
presented  in  Table  3.4.  The  effectiveness  of  these  checks  is 
demonstrated  with  Table  3*5  by  observing  K1-K4,  before  and  after  the 
checks  are  implemented. 

SOLVING  FOR  THE  RADIAL  FIELD.  The  radial  electric  field  is  solved 
by  applying  a  fourth  order  Runge-Kutta  technique  to  the  differential 
equation  for  the  radial  electric  field  (Eq  3.63).  The  independent 
variable  is  retarded  time.  Equation  3.66  becomes 

dE^ 

—  =  f(T,E^)  with  r=const  (3.75) 

dx 


K1 

K2 

K3 

K4 

J/SIGMA 

E-FIELD 

(Volts/Meter) 

9524. 

9394. 

28.7 

25.6 

29.9 

11.8 

9524. 

9418. 

25.7 

16.5 

31.5 

-20.5 

9524. 

9435. 

31  .k 

0.61 

59.6 

-155. 

9524. 

9435. 

92.2 

-81.4 

299. 

-1245. 

9524. 

93’5. 

707. 

-969. 

3126. 

-1.5E4 

9525. 

7622. 

9^*76. 

-1 .5E4 

4.9E4 

-2.5E5 

9525. 

-2.2E4 

1.7E5 

-2.8E5 

9.3E5 

-4.8E6 

9527. 

-5.8E5 

Table 

3,4  Runge- 

-Kutta  Without 

Automatic 

Step  Reduction 

K1 

K2 

K3 

K4 

J/SIGMA 

E-FIELD 

(Volts/Meter) 

9524  . 

9394. 

28.7 

25.6 

29.9 

1 1 .8 

9524. 

9418. 

10.6 

8.97 

10.3 

7.08 

9524. 

9440. 

7.60 

6.52 

7.60 

4.90 

9524. 

9455. 

5.67 

4.69 

5.80 

3.00 

9524  . 

9465. 

4.35 

3.21 

4.63 

1.09 

9525. 

9472. 

1.13 

1.10 

1.12 

1.07 

9525. 

9478. 

0.54 

0.53 

0.54 

0.52 

9526. 

9480 

Table  3 

.5  Runge-Kutta  With 

Automatic 

Step  Reduction 

5i 
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The  radial  electric  field  is  solved  by  stepping  along  retarded  time  with 
position  held  constant.  The  ray  below  symbolizes  this  displacement. 


where  '^n''^n+1  Compton  current  for 

®n'°n+1  -conductivity  for  positions  ’'^^'^n+l 
-radial  electric  field  for  position 
-unknown  radial  electric  field  for  position 

EMPFRE  generates  values  of  the  current  and  conductivity  for 
positions  and  but  the  Runge-Kutta  method  requires  an  additional 

value  at  the  mid  point.  The  current  and  conductivity  are  assumed  to  be 
linear  across  At;  therefore,  the  mid  point  values  are  the  linear 
averages  of  the  values  for  positions  and  If  the  retarded  time 

step  is  reduced,  additional  values  are  generated.  Equations  that 
generate  the  estimates  of  the  change  in  E^,  over  the  interval  become 


K1=-  — -(Vo^-  E^) 
^o 

At 

K2=-  _  . _  + 

"o  L  2 

At  RJn+Jn^l) 

K3  =  -  —•  -  + 

^o  L  2 


•  (Ejj+K1/2) 


■  (E^+K2/2) 


(3.76) 


(3.77) 


(3.78) 
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(3.79) 


At 

K4  =  - 

"o 

The  unknown  radial  electric  field  at  is  approximated  by 

^n+1  =En+[K1+2  (K2+K3)+KiJ]/6  (3.80) 

The  radial  electric  field  is  assumed  to  be  zero  for  all  R'  at  t=0. 

SOLVING  FOR  THE  TRANSVERSE  ELECTRIC  FIELD.  The  transverse  electric 
fields  are  solved  by  applying  a  fourth  order  Runge-Kutta  technique  to 
the  differential  equation  for  the  theta  or  phi  electric  fields  (Eqs  3.6^ 
or  3.65).  The  theta  and  phi  electric  fields  are  considered  together  in 
this  section  because  the  differential  equations  are  identical  in  form 
and  the  method  outlined  in  this  section  applies  to  either. 

The  independent  'ariable  is  the  spatial  position  R.  Equation  3.66 
becomes 


6Ey 

-  =  f{r,E^)  (3.81) 

6r 
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A  transverse  electric  field  is  solved  by  stepping  along  the  slant 


range  with  retarded  time  held  constant.  The  ray  below  symbolizes  this 
displacement : 


T 


■"n 

^n+1 

°n'‘^n*^n 

°n+1**^n+1'^n+1 

r 

where 

JntJn+1  -a  transverse  Compton  current  for  r^,  r^^^  at  t’ 

On'^’n+I  -conductivity  for  r^,  r^j^^  at  t* 

Efj  -a  transverse  electric  field  for  r^ 

-unknown  transverse  electric  field  for  r^^^^ 

EMPFRE  generates  values  of  the  current  and  conductivity  for 
positions  r^^  and  but  the  Runge-Kutta  method  requires  an  additional 

value  at  the  mid  point.  The  current  and  conductivity  are  assumed  to  be 
linear  across  Ar;  therefore,  the  midpoint  values  are  the  linear  averages 
of  the  values  for  the  positions  r^  and  If  the  spatial  step  is 

reduced,  additional  values  are  generated.  Equations  that  generate  the 
estimates  of  the  change  in  E^  over  the  integral  Ar  become 


K1=  -Ar*  _  + 
r 


1  Mo*c*o„ 


Uo’C 


■V- 

2  2 


(3.82) 


K2=-Ar-. 


1  U  •  C*  ( 0  •♦•0  4  )  u  •  c 

'  ^o  '  n  n+1'  ,  .0/7.1  )  (0  Q-3) 

_  +  _ , _  •  (E„+K1/2)+ _ —^'^n^'^n+V 


r^+Ar/2 
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K3=-Ar 


r^+Ar/2 


_ •  (V^2/2)+ _ (Jn-"Jn+l)  (3.8i<) 

H  4 


r 


n+1 


U  •  C*  0 

n+ 


•Ie^.K3U - 

2  2 


(3.85) 


The  unknown  transverse  electric  field  for  is  approximated  by 


En^1=En+(K1+2- (K2+K3)+K4)/6 


(3.86) 


The  transverse  electric  field  is  assumed  to  be  zero  above  the  absorption 
region  for  all  t. 
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COMPUTATIONAL  PROCEDURE.  The  electron  mobility  is  a  function  of 
the  total  electric  field.  The  conductivity  is  a  function  of  electron 
mobility;  therefore,  the  three  differential  equations  for  the  electric 
field  are  dependent  and  must  be  solved  simultaneously.  The  procedure 
for  solving  for  the  electric  field  is  symbolized  in  Figure  3.3-  An 
estimate  of  the  electric  field  at  ^n+1  made  to  calculate 

the  conductivity.  The  conductivity  at  is  based  on  the  total 

electric  field  (E^)  from  r^^^,  t^.  This  direction  was  chosen  because 
the  number  of  retarded  time  steps  is  large  compared  to  the  number  of 
positional  steps. 

The  fields  along  are  solved  for  all  r.  The  fields  along 
have  been  solved  from  the  top  of  the  absorption  region  down  to  r^.  The 
transverse  fields  at  r^^^,  are  solved  from  the  fields  at  r^^, 

The  radial  field  at  are  solved  from  the  radial  field  at 

T|^.  This  procedure  is  repeated  until  the  lower  limit  of  the  absorption 
region  is  reached.  The  transverse  components  of  the  EMP  will  attenuate 
as  1/r  as  it  propagates  to  the  ground. 


59 


FOURIER  TRANSFORM 


The  Fourier  transform  subroutine  solves  the  following  integrals 


using  Simpson’s  rule: 


E(t)  E(t) 

g(a)=  /  — — .  •COS(oiT)dT-i  /  -  •SIN(aT)dT 

-00  v  2ir  -00  /  2Tt 


(3.89) 


The  magnitude  of  the  Fourier  transform  is  found  by  multiplying  by 
its  complex  conjugate  and  taking  the  square  root  of  the  product: 


|g(a)|  =  [g(a)*g  (a) 


(3.90) 


A  relationship  between  the  energy  associated  with  a  frequency  band 
and  the  Fourier  transform  is  found  by  applying  Parseval's  theorem: 


/  g(a)]  ^d  a=  /  |e(t)|  ^di 


(3.91 ) 
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The  energy  (e)  associated  with  the  electromagnetic  pulse  is  equal  to 
”  |e(t)| 2 

e  =  /  -  dT  (3.92) 

o  n 

where  n  = 

Substituting  this  expression  of  energy  into  Parseval's  theorem 
relates  the  Fourier  transform  to  energy: 

“  jg(a)|2 

e  =  /  -  da  (3.93) 

o  n 

The  continuous  Fourier  transform  is  evaluated  at  30  discrete  values 
of  angular  frequency.  The  energy  associated  with  a  frequency  is  found 
by 

|s(a)p 

Energy  (J-sec/m^)  =  _  (3.9^) 

n 

To  save  execution  time,  the  late  time  electric  fields  are  not 
calculated.  Terminating  the  integral  for  the  Fourier  transform  at  the 
last  value  for  the  electric  field  E(Tjg),  would  produce  an  artificial 
transient.  The  abrupt  termination  would  effect  the  frequency  content. 
The  Fourier  transform  subroutine  fits  a  sine  function  for  the  electric 
field  from  E(Tj^)  to  an  artificial  termination  at  E(tj,)=0.  The  slope  of 


the  sine  function  is  matched  to  the  slope  the  electric  field  at  The 
angular  frequency  for  the  sine  function  (oi)  is  found  by 


where  At  -retarded  time  step 

A  -amplitude  of  the  sine  function 


For  T  >  T^,  the  electric  field  is  set  to 
E(t)  =  A[1-SIN(aJ*(T-T^))]  (3-96) 

The  upper  limit  of  the  integral  becomes 

T  =  T  +  it/(2u))  (3*97) 

where  ECt^)  =  0  and  dE(T) 

dX  Tg 
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IV.  ANALYSIS  OF  NUMERICAL  TECHNIQUE 


The  solution  of  Maxwell's  equations  for  high  altitude  EMP  in  EMPFRE 
involves  several  trade-offs  between  accuracy  and  computer  resources. 
The  purpose  of  EMPFRE  is  to  show  the  effects  burst,  parameters  such  as 
geometry  and  gamma  output,  on  the  frequency  content  of  the  EMP.  A  small 
accumulative  error  is  accepted  because  the  basic  cause  and  effect 
relationships  will  be  insensitive  to  the  accumulative  error  associated 
with  EMPFRE.  An  effort  is  made  to  identify  exceptions  to  this  rule. 
Parametric  studies  are  made  to  determine  optimum  compromises  for  each 
trade-off.  A  review  of  these  studies  is  made  here  to  provide  insight 
about  the  error  accepted  in  EMPFRE.  The  input  parameters  for  the  base 
case  are  listed  below; 


X=0  Meters 

(^.  1 ) 

Y=0  Meters 

(‘*.2) 

Z=0  Meters 

(4.3) 

HOB=100  Kilometers 

(4.4) 

Y.^=0.1  Kilotons 

(4.5) 

E^=2  MeV 

(4.6) 

8=5.5x10“^  Webers/Meter 

(4.7) 

Dip  Angle=0  Degrees 

(4.8) 

NDELR=50 

(4.9) 

DELT=5x10"’°  Seconds 

(4.10) 

These  parameters  are  held  constant  while  the  parameter  in  question  is 
varied.  The  acculative  error  is  estimated  by  using  a  sample  HEMPS  (High 
Altitude  EMP  Program  used  at  the  Air  Force  Weapons  Lab)  solution  as  a 


test  case.  The  comparison  to  HEMPB  is  made  in  Appendix  A. 

COMPTON  ELECTRON  RANGE 

The  Compton  electron  range  determines  the  secondary  electron 
production  rate  in  EMPFRE.  The  program  assumes  a  linear  energy  loss 
over  the  life  of  the  Compton  electron.  Therefore,  the  secondary 
electron  production  rate  is  proportional  to  the  total  local  Compton 
current.  With  this  model,  changing  the  effective  electron  range  for  a 
fixed  Compton  electron  energy  will  impact  the  secondary  electron 
production  rate.  The  physical  effect  being  observed  by  varying  the 
Compton  electron  range  is  changing  the  secondary  electron  production 
rate . 

The  mean  forward  range  for  the  Compton  electrons  is  approximated  by 
a  relationship  from  the  Radiological  Health  Handboox  ^Ref  16:29).  The 
estimate  of  the  electron  range  for  a  2  MeV  gamma  ray  is  4.06  meters  at 
sea  level.  The  Compton  electron  range  is  varied  from  2.03  to  6.09 
meters.  Figure  <4.1  shows  the  variation  in  the  resulting  EMP.  Figure 
4.2  shows  the  variation  introduced  into  the  frequency  content. 

The  estimate  of  the  mean  forward  range  of  the  Compton  electron 
impacts  the  secondary  electron  density  and  the  conductivity.  The 
lifetime  of  the  Compton  electron  is  calculated  from  the  initial  electron 
velocity  and  the  forward  range;  therefore,  an  increase  in  the  forward 
range  will  increase  the  Compton  electron  lifetime.  The  number  of 
secondary  electrons  produced  per  Compton  electron  is  fixed  at  l/33eV. 
The  secondary  electron  production  rate  is  inversely  proportional  to  the 
Compton  electron  life-time.  An  increase  in  the  Compton  electron  life¬ 


time  will  decrease  the  secondary  electron  production  rate  and  density  at 


Figure  4.1  Electric  Field,  Compton  Electron  Pange 
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early  times.  The  resulting  decrease  in  conductivity  will  cause  the  EMP 
to  rise  faster  and  decay  slower.  The  resulting  electric  fields  are 
larger  and  the  pulse  is  broadened. 

The  change  in  the  electric  fields  has  an  impact  on  the  frequency 
profile  of  the  EMP.  The  threshold  frequency  is  defined  as  the  frequency 
where  the  energy  begins  to  drop  from  its  constant  value  for  low 
frequencies.  The  energy  increases  with  Compton  range  but  the  threshold 
frequency  decreases  due  to  broadening  of  the  pulse.  An  error  in  the 
Compton  electron  range  has  a  large  impact  on  the  resulting  EMP. 

SPATIAL  AND  TIME  INCREMENTS 

A  fourth  order  Runge-Kutta  numerical  technique  is  used  to  solve  the 
differential  equations  for  the  electric  fields.  The  application  of  this 
method  is  discussed  in  Chapter  III.  The  solution  beco'aes  unstable  when 
the  spatial  (for  the  transverse  field)  or  time  (for  the  radial  field) 
increments  become  large.  This  problem  is  corrected  by  internal  checks 
that  halve  the  increment  if  instabi lities  are  predicted.  The  remaining 
error  is  due  to  the  discrete  parameter  updates. 

The  key  parameter  is  the  electron  mobility.  The  electron  mobility 
is  a  function  of  the  electric  field  and  altitude.  The  electric  field  is 
approximated  by  the  electric  field  from  the  last  time  step  at  that 
position.  See  Figure  3.6  for  a  stepwise  procedure  in  the  solution  of 
the  electric  field.  The  variation  in  the  EMP  associated  with  increasing 
the  time  increment  is  demonstrated  by  Figure  i4.3.  The  time  increment  is 
varied  between  2.5x10"^^  and  2x^0'''^  seconds.  Figure  shows  the 
associated  variation  in  the  frequency  content. 
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Figure  4.4  Frequency  Profile,  Time  Step 


Increases  in  the  spatial  increment  produce  an  error  because 
parameters  such  as  Compton  currents,  secondary  electron  density  and 
electron  mobility  will  vary  across  the  spatial  increment.  The  variation 
in  the  EMP  associated  with  increasing  the  spatial  increment  is 
demonstrated  by  Figure  4.5.  The  spatial  increment  is  varied  between  25 
and  100  steps.  Figure  4.6  shows  the  associated  variation  in  the 
frequency  profile. 


LORENTZ  MODEL 

The  Lorentz  model  assumes  the  secondary  electrons  are  born  at 
thermal  equilibrium  with  the  electrons  drifting  in  the  local  electric 
field.  This  approximation  is  imposed  to  simplify  the  calculations  but 
causes  EMPFRE  to  underestimate  the  rise  time  and  peak  of  the  EM^*.  As  a 
result,  the  energy  associated  with  the  high  fr'equencies  are 
underestimated. 

The  secondary  electrons  are  born  at  about  10  eV  (Ref  2)  and  take 

q 

more  than  1x10  seconds  to  therraalize  from  this  energy  at  an  altitude 
of  20  kilometers.  By  assuming  the  velocities  are  isotropic,  the 
secondary  electrons  are  effectively  at  a  higher  temperature  than 
proposed  by  the  Lorentz  model.  A  more  detailed  discussion  of  this  issue 
is  presented  in  Chapter  VI.  The  lower  temperature  predicted  by  the 
Lorentz  model  over-estimates  the  electron  mobility  because  the  collision 
cross  sections  vary  with  electron  energy.  In  the  Lorentz  model,  the 
secondary  electron  mobility  ranges  from  2.3'^  to  0.195  ra^/Csec-voltl  at 
sea  level.  A  crude  way  to  demonstrate  the  variation  in  the  EMP 
associated  with  the  depressed  secondary  electron  temperature  is  to  set 
an  upper  limit  on  electron  mobility.  The  base  case  uses  the  Lorentz 


Figure  4.5  Electric  Field,  Spatial  Step 
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Figure  4.6  Frequency  Profile,  Spatial  Step 


Model  and  the  maximum  mobility  at  sea  level  is  varied  from  the  Lorentz 
model  to  a  maximum  mobility  of  0.02  m^ /(sec- vo  1 1)  and  adjusted  for 
altitude.  The  variation  in  the  EMP  associated  with  the  increase  in 
secondary  electron  temperature  is  demonstrated  in  Figure  ii.7.  Figure 
4.8  shows  the  associated  variation  in  the  frequency  content. 

TRUNCATION  OF  ELECTRIC  FIELD 

EMPFRE  calculates  the  electric  field  for  several  shakes  and  then 
terminates  the  calculation.  Test  runs  indicate  the  time  step  size  could 
be  increased  after  the  first  few  shakes  without  introducing  excessive 
error.  This  approach  would  be  appropriate  for  examining  the  electric 
field  but  not  for  studying  the  high  frequency  characteristics.  .A 
decrease  in  the  sampling  rate  would  cause  aliasing  in  an  unpredictable 
manner.  Termination  of  the  pulse  is  accomipl  ished  by  merging  a  sine 
function  to  it  as  discussed  in  Chapter  III. 

The  angular  frequency  of  the  sine  function  is  chosen  to  match  the 
slope  of  the  pulse  at  termination.  The  impact  of  this  technique  on  the 
Fourier  transform  is  discussed  in  Appendix  A, 
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Figure  4.8  Frequency  Profile,  Mobility  Limit 


V.  PARAMETEHIC  STUDIES 


The  gamma  ray  pulse  from  a  device  rises  exponentially.  The 
parameter  for  characterizing  the  rise  will  vary  due  to  changes  in 
compression  and  boosting.  The  details  of  the  gamma  ray  pulse  history 
are  not  modeled  due  to  classification.  EMPFRE  approximates  the  gamma 
ray  pulse  history  with  a  single  parameter  for  exponential  rise  and  a 
second  parameter  for  decay. 


(t)  =  c 


• 

.  1  + 


EXP (ax) 


EXP(B* (t-TP)) 


(5.1) 


where  a  -parameter  to  characterize  exponential  rise 
3  -parameter  to  characterize  exponential  decay 
TP  -time  of  peak 
C  -yield  normalization  constant 


For  retarded  time  less  than  TP,  the  gamma  output  of  the  device  rises 
exponentially . 


Y.^(x)  =  C*EXP(aT)  (5.2) 

For  retarded  time  greater  than  TP,  the  gamma  output  of  tne  device  decays 
exponentially. 


Y^(t)  =  K*EXP((a-3)-T)  (5-3^ 

This  pulse  shape  is  used  for  the  parametric  studies.  'Jnles.:  specLfied, 
the  parameters  for  the  burst  are; 

a  rl .0E9/3econd  (5.9) 


P  = 1 . 5E9/3econd 


(5.5' 


TP  =1.0E-8  second 


(5.6) 


E  ^=1.5  MeV  (5.7) 

Yy  =0.1  kiloton  (5.8) 

HOB  =  100  kilometers  (5.9) 

BANGLE  =0  degrees  (5.10) 

TARGET  LOCATION: 

X  =0  meters  (5.11) 

Y  =0  meters  (5.12) 

Z  =0  meters  (5. tj  ) 

GAMMA  PULSE  HISTORY 


The  alpha  (a)  of  the  gamma  pulse  impacts  the  EMP  rise  time  and  peak 
field.  The  gamma  ray  rise  time  determines  the  Compton  current  history 
which  directly  impacts  the  EMP.  Another  mechanism  involves  the  phasing 
between  the  Compton  current  and  conductivity.  The  EMP  rises  with  the 
Compton  currents  until  the  air  conductivity  becomes  large.  The  Coraptoh 
electrons  have  a  time  of  flight  greater  than  10  shakes;  therefore,  the 
Compton  currents  continue  to  rise  while  the  gamma  ray  pulse  decays.  The 
decay  of  the  EMP  is  due  to  the  increase  in  conductivity.  The  rise 
characteristics  of  the  gamma  pu’se  determines  how  far  ttie  electric  field 
rises  before  being  choked  off  by  the  build-up  of  conductivity. 

The  impact  of  a  change  in  the  parameter  alpha  [Eq.  5.1]  on  the 
electric  field  is  demonstrated  by  Figures  5.1,  5-3,  5.5  and  5.7.  Thi' 
impact  on  the  spectral  characteristics  are  indicated  by  Figures  5.’, 
5.^,  5.6  and  5.8.  The  examples  are  grouped  with  equal  values  of  TP. 
The  TP  is  chosen  to  reduce  the  effects  of  the  function  starting  at  a 


non-zero  value. 
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Figure  5.5  Electric  Field,  TP=1.0  Shake,  Beta=1 .5  Alpha 
Alpha=5.0x10°,  1.0x10^  and  2.5x10^  (sec"^ 
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GAMMA  RAY  ENERGY 


The  gamma  ray  spectrum  from  a  nuclear  device  will  be  a  function  of 
the  mix  of  nuclear  reactions  and  the  composition  of  the  device.  The 
purpose  is  to  demonstrate  the  relationship  between  mean  gamma  ray  energy 
and  the  spectral  characteristics  of  the  EMP.  No  attempt  is  made  to 
model  the  EMP  from  a  spectrum  of  gamma  rays.  This  would  require  a 
multi-group  gamma  ray  calculation  or  a  model  to  account  for  spectrum 
hardening.  Spectrum  hardening  is  the  effect  of  having  the  mean  energy 
of  the  gamma  spectrum  increase  as  the  low  energy  gamma  rays  are  filtered 
out  faster  than  high  energy  gamma  rays. 

Changing  the  gamma  ray  energy  changes  the  Compton  electron  density 
profile  as  discussed  in  Chapter  II.  Increasing  the  energy  of  the  gamma 
rays  causes  the  peak  of  the  deposition  region  to  shift  to  lower 
altitudes.  The  electron  mobility  is  inversely  proportional  to  pressure; 
therefore,  the  electron  mobility  is  reduced  at  the  peak  of  the 
deposition  region.  The  resulting  reduction  in  conductivity  lead  to 
larger  electric  fields. 

The  Compton  electron  dynamics  are  effected  by  the  gamma  ray  energy. 
An  increase  in  the  gamma  ray  energy  will  cause  increases  in  the  Compton 
electron  energy,  initial  velocity  and  range.  The  Compton  electrons  will 
generate  secondary  electrons  faster  which  enhances  the  build-up  of  the 
secondary  electron  density.  An  increase  in  the  gamma  energy  will 
decrease  the  mean  deflection  angle  of  the  Compton  electron  (Table  3-?^* 
This  along  with  the  higher  intial  velocity  increases  the  Compton 
current.  The  net  effect  from  Compton  electron  dynamics  leads  to  faster 
rise  times  and  higher  peak  electric  fields. 


The  yield  is  set  at  0.1  kilotons.  This  has  the  effect  of  reducting 
the  number  of  gamma  rays  emitted  from  an  increase  in  gamma  ray  energy. 
This  will  reduce  the  Compton  currents.  The  impact  from  changing  the 
energy  of  the  gamma  ray  with  this  approach  is  demonstrated  by  Figure 
5.9.  The  variation  in  the  spectral  characteristi cs  is  indicated  in 
Figure  5. 10. 

Another  way  to  see  the  effects  of  changing  the  mean  energy  is  to 
fix  the  number  of  gamma  rays  emitted.  The  impact  from  changing  the 
energy  with  this  approach  is  demonstrated  by  Figure  11.  The  10  MeV 
gamma  ray  pulse  is  largest  but  close  to  the  5  MeV  pulse.  The  effects 
due  to  the  Compton  electron  dynamics  are  being  cancelled  by  gamma 
absorption  from  pair  production.  The  number  of  gamma  rays  is  based  on 
0.1  kiloton  yield  of  1.5  MeV  gamma  rays.  The  variation  in  the  spectral 
characteristics  is  indicated  in  Figure  5.12. 
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Figure  5«10  Frequency  Profile,  Gamoa  Ray  Energy 
YleldaO.1  Kllotons 
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Figure  5.12  Frequency  Profile,  Gamma  Ray  Energy 
YieldsE  (MeV)- 0.0667  Kilotons 


YIELD 


The  yield  impacts  the  Compton  currents  and  the  conductivity 
histories.  Raising  the  yield  by  an  order  of  magnitude  (0.1  Kt  1.0  Kt) 
produced  an  increase  of  only  25%  in  the  peak  electric  field.  This  is 
due  to  the  effects  of  saturation.  The  concept  of  saturation  can  be 
understood  by  analyzing  the  differential  equation  for  the  transverse 
electric  fields. 

6  E.J.  E.J1  u*  C 

-  =  -  —  -  - '(J^+a^E^)  (5.1^) 

6r  r  2 

For  electric  fields  greater  than  the  saturation  field  the  conduction 
current  (a*E)  would  be  larger  than  the  Compton  current.  Ignoring  the 
radial  divergence,  the  electric  field  will  remain  less  than  the 
saturation  field  (-J/a).  The  air  conductivity  is  increasing  with  the 
Compton  currents  to  hold  the  peak  field  down.  The  increase  in  the 
electric  field  is  due  to  a  change  in  the  saturation  field. 

At  lower  yields  the  electric  field  is  not  saturated  and  the 
conductivity  remains  low.  The  result  Is  a  lower  peak  electric  field  but 
higher  electric  fields  for  time  greater  than  3  shakes. 

Increasing  the  yield  to  saturation  leads  to  a  shorter  pulse  with  a 
faster  rise  and  decay.  This  Increases  the  energy  associated  with  the 
high  frequencies  with  little  change  in  the  low  frequencies.  A  further 
increase  in  yield  (to  =  1.0  kiloton)  increased  the  low  frequencies 
with  little  change  in  the  high  frequencies.  The  impact  of  a  change  in 
yield  on  the  electric  field  is  demonstrated  by  Figure  5.13.  The  impact 
on  the  frequency  profiles  is  in  Figure  5. 14. 
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Figure  5.14  Frequency  Profile,  GamDa  Ray  Yield 


HEIGHT  OF  BURST 


Changing  the  height  of  burst  has  an  effect  similar  to  changing  the 
yield.  The  Compton  current  and  conductivity  are  decreased  for  larger 
heights  of  burst.  The  deposition  profiles  for  various  HOBs  are  given  in 
Chapter  II.  At  low  HOB  (75  Kra  and  100  Km),  the  fields  are  limited  by 
the  saturation  voltage;  therefore,  the  peak  voltage  is  relatively 
insensitive  to  decreasing  HOB.  For  large  HOB  (200  Km),  the  conductivity 
remains  low  and  the  electric  field  doesn't  saturate.  As  a  result,  the 
electric  fields  are  larger  for  time  greater  than  3  shakes.  Figure  3.15 
demonstrates  the  Impact  on  the  electric  fields.  The  frequency  profiles 


DIP  ANGLE 


The  dip  angle  of  the  geomagnetic  field  will  change  the  trajectory 
of  the  Compton  electrons.  The  phi  component  of  the  Compton  current  is 
largest  for  the  dip  angle  equal  to  zero.  The  phi  component  of  the 
Compton  current  will  go  to  zero  as  the  dip  angle  goest  to  90^^.  The 
theta  component  is  zero  when  the  dip  angle  equals  zero.  The  theta 
component  reaches  a  maximum  when  the  dip  angle  is  equal  to  ■45'^.  Then 
the  theta  component  goes  to  zero  as  the  dip  angle  goes  to  90°. 

When  the  dip  angle  is  at  90°  (parallel  to  the  gamma  ray 
trajectory),  there  are  no  transverse  currents  generated.  In  this  case, 
the  Compton  currents  do  not  form  a  radiating  magnetic  dipole.  The 
resulting  EMP  would  be  due  to  the  electric  dipole  formed  by  the  radial 
currents.  The  electric  fields  radiating  from  the  electric  dipole  are 
small  and  are  lower  in  frequency  compared  to  those  from  the  magnetic 
dipole  for  the  case  of  high  altitude  EMP.  The  phi  component  of  the 
Compton  current  makes  the  largest  contribution  to  the  magnetic  dipole. 
The  radiated  electric  fields  go  to  zero  as  the  dip  angle  goest  to  90°. 

The  components  of  the  Compton  current  are  a  function  of  dip  angle 
but  the  conductivity  profile  will  remain  the  same.  As  a  result,  the 
peak  electric  field  does  not  shift  with  dip  angle.  The  effect  of  dip 
angle  on  the  electric  field  is  demonstrated  by  Figure  5.17.  The 
frequency  profiles  are  in  Figure  5,18. 


98 


VI.  CONCLUSIONS  AND  RECOMMENDATIONS 


CONCLUSIONS 

The  EMP  from  high  altitude  nuclear  detonations  is  modeled  by 
EMPFRE.  Constraints  on  the  project  made  it  impractical  to  refine  the 
model  or  reduce  the  error  by  imposing  finer  steps  for  computations. 
These  limitations  do  not  invalidate  the  results  of  the  parametric 
studies  performed.  This  project  documents  basic  relationships  between 
burst  parameters  and  the  resulting  EMP.  "’he  following  conclusions  are 
drawn  for  the  parametric  studies. 

1)  The  energy  in  the  EMP  is  enhanced  by  a  large  alpha.  This  is 
due  to  the  phasing  between  the  Compton  current  and  conductivity.  An 
increase  in  alpha  allows  the  Compton  currents  to  get  larger  before  the 
conductivity  gets  large.  The  result  is  faster  rise  time  and  larger  pea’K 
fields.  The  energy  associated  with  100  MHs  is  several  orders  of 
magnitude  lower  for  alpha= 1.5X10^/3econd  than  for  alpha=2.5X10^/second. 

?.)  Tne  energy  in  the  EMP  is  enhanced  by  higher  gamma  ray  energies; 
however,  a  limit  is  reached  where  pair  production  counteracts  the 
benefits  due  to  Compton  electron  dynamics.  The  energy  in  the  EMP  is 
increased  by  more  than  an  order  of  magnitude  when  the  gamma  ray  energy 
is  changed  from  1  MeV  to  2,5  MeV.  Increases  in  the  mean  energy  above  5 
MeV  results  in  a  reduction  of  the  energy. 

3)  The  energy  in  the  EMP  is  enhanced  by  an  increase  in  yield. 
This  is  particularly  true  for  frequencies  above  10  MHz.  The  increase  in 
yield  produces  faster  rise  times  and  decay  times.  Higher  yields  produce 
larger  peak  fields  and  narrower  pulses.  Increases  in  yield  passed  that 
needed  to  saturate  the  electric  field  produce  limited  changes  in  the 
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energy  in  the  EMP. 

^4)  Decreasing  the  angle  of  intersection  between  the  geomagnetic 
field  and  the  slant  range  decreases  the  energy  associated  with  the 
electric  field.  Varying  this  angle  between  90°  and  60°  has  little 
effect  on  the  resulting  EMP.  Varying  this  angle  between  60°  and  30° 
results  in  a  large  reduction  in  the  EMP.  Tnis  effect  is  due  to  a 
reduction  of  the  phi  component  of  the  Compton  current. 

RECOMMENDATIONS 

SPATIAL  STEP  REDUCTION.  The  behavior  of  the  late  time  electric 
fields  could  be  improved  by  increasing  the  number  of  spatial  steps  from 
5U  to  at  least  100.  Figure  9.5  shows  the  improvement  made  for  the  bas.' 
case.  A  larger  number  of  steps  may  be  required  for  some  cases.  Thi: 
effects  on  the  frequency  profile  were  small  for  the  base  case  but 
confidence  would  be  increased  in  other  results  if  this  recommendation  is 
implemented . 

REFINEMENTS  TO  LORENTZ  MODEL.  The  Lorent:i  model  is  an 

approximation  used  to  model  secondary  electron  behavior.  The  secondary 
electrons  are  assumed  to  be  created  at  the  equilibrium  temperature  of 
the  electrons  in  the  local  electric  field.  The  potential  impact  of  this 
approximation  is  demonstrated  by  Figure  9.7. 

At  early  times,  the  growth  in  the  Compton  current  density  is 
matched  by  a  growth  in  the  e lectroraagnet io  fields  generated.  The 
simultaneous  growth  of  secondary  electron  density  increases  air 
conductivity  and  limits  the  growth  of  the  electromagnetic  field.  The 
peak  field  strength  and  frequency  profile  are  sensitive  to  the  phasing 
of  the  Compton  current  density  and  conductivity  growth. 


Secondary  electrons  are  born  with  temperatures  above  the 
equilibrium  temperature.  For  dry  air,  the  electron  mobility  is  reduced 
by  an  increase  in  the  temperature.  A  derivation  of  the  relationship 
between  electron  temperature  and  mobility  is  presented  by  C.  E.  Baum 
(Ref  1).  The  Lorentz  model  overestimates  the  secondary  electron 
mobility  from  birth  to  thermalization. 

Two  approaches  are  identified.  The  first  approach  is  discussed  by 
L.  W.  Seiler  (Ref  15:2^).  The  secondary  electron  collision  frequency  is 
assumed  to  be  linear  function  of  time  between  birth  and  thermalization. 
Another  approach  is  used  by  Karzas  and  Latter  (Ref  10).  In  the 
reference.  Age  Theory  is  used  to  find  the  average  electron  mobility  for 
a  given  alpha.  This  is  the  same  approach  used  to  find  single  group 
cross  sections  for  systems  of  neutrons.  The  last  approach  appears  to  be 
more  flexible  for  application  in  EMPFRE. 

SPECTRUM  HARDENING.  Incorporating  a  model  to  treat  spectrum 
hardening  is  proposed  as  a  follow-up  thesis.  The  first  step  is  the 
development  of  a  method  to  determine  an  effective  gamma  ray  energy  for  a 
given  gamma  ray  spectrum.  Important  cha rac ter is ti cs  to  consider  are 
MFP,  Compton  electron  energy  and  Compton  electron  deflection  angle.  The 
second  step  is  to  model  the  shift  in  effective  gamma  energy  as  the  pulse 
propagates  through  the  atmosphere.  Lt  A.  Dooley  (GNE-3'1M)  worked  on  a 
thesis  to  model  the  impact  on  gamma  ray  spectral  characteristics  due  to 
attenuation  by  the  atmosphere.  The  results  from  this  work  may  provide 
the  data  needed  to  develop  a  model  for  spectrum  hardening. 


103 


BIBLIOGRAPHY 


1.  Baum,  C.  E.  "The  Calculation  of  Conduction  Electron  Parameters  in 

Ionized  Air,"  Electromagnetic  Pulse  Theoretical  Notes:  Note  6, 
Kirtland  AFB,  New  Mexico:  Air  Force  Weapons  Laboratory,  (Apcfl 
1967) 

2.  ,  "Electron  Thermalization  and  Mobility  in  Air,"  Electromagnetic 
Pulse  Theoretical  Notes:  Note  12,  Xirtland  AFB,  New  Mexico:  Air 
Force  Weapons  Laboratory,  (April  1967) 

3.  Broad,  W.  J.  "Nuclear  Pulse  (I):  Awakening  to  the  Chaos  Factor," 

Science:  1009-1012  (May  1931) 

4.  ,  "Nuclear  Pulse  (II):  Ensuring  Delivery  of  the  Doomsday 

Signal,"  Science :  1116-1120  (June  1931) 

5.  Erkkila,  J.  H.  "Calculations  of  the  EMP  From  High  Altitude  Nuclear 
Detonations,"  Electromagnetic  Pulse  Theoretical  Notes :  Note  26. 
Kirtland  AFB,  New  Mexico:  Air  Force  Weapons  Laboratory,  rAp'’il 
1967  ) 

6.  ,  "Prompt  Gamma-Ray  Fluxes  and  Energy  Deposition  in  an 

Exponential  Atmosphere,"  Electromagnetic  Pulse  Theoretical  Notes : 
Note  17,  Kirtland  AFB,  New  Mexico:  Air  Force  Weapons  Laboratory, 
t'Aprll  1967  ) 


7.  ,  Lecture  Materials  distributed  in  NE-625.  School  of 

Engineering,  Air  Force  Institute  of  Technology  (AU),  Wright- 
Patterson  AFB  OH,  (1983) 

3.  Evans,  R.  D.  The  Atomic  Nucleus.  New  York;  McGraw  Hill  Book  Co. 
(1966) 

9.  Karza.s,  W.  J.  and  Richard  Ufitter.  Detection  of  the 
Electromagnetic  Radiation  From  Nuclear  Explosions  in  Space," 
Physical  Review,  Vol  137,  Number  5B  (March  1965) 


10.  ,  Air  Conductivity  Produced  by  Nuclear  Explosions.  MEMORANDUM 

RM-3671-PR,  RAND  Corporation,  'May  1963'' 

11.  Longley,  H.  J.  and  Conrad  L.  Longmire.  E 1 ectron  Mobl 1 i ty  and 

Attachment  Rate  in  Moist  Air.  MHC-M-222,  Mi.ssion  Research 
Corporation,  Santa  Barbara,  CA  '’December  1975) 

12.  Longmire,  Conrad  L.  "On  the  E 1  ectroiaagne  t  Lc  Pulse  Produced  by 

Nuclear  Explosions,"  IEEE  Transactions  on  Antennas  and  Propagation, 
AP-26:  3-13  (January  19781 


1 


13.  Milne,  W.  E.  Numerical  Solutions  of  Differential  Equations.  New 

York:  Dover  Publications  Inc.  (1970) 

14.  Pettus,  E.  and  W.  F.  Crevier.  Analytic  Representation  of  Electron 

Mobility  and  Attachment  Data  In  Dry  and  Moist  Air  from  van  Lint’ s  HIFX 
Experiments.  MRC-R-ST^  ^anta  Barbara,  CA:  Mission  Research 
Corporation  (July  1980) 


15.  Seiler,  L.  W.  A  Ca leu lationa  1  Model  for  High  Altitude  EMP.  MS 

thesis.  School  of  Engineering,  Air  Force  Institute  of  Technology  (AU), 
Wright-Patterson  AFB  OH  (March  1975) 

16.  U.  S.  Department  of  Health  Education  and  Welfare.  Radiological 

Health  Handbook.  Revised  edition  (January  1970) 


APPENDIX  A 


VERIFICATION  ^  EMPFRE 

FOURIER  TRANSFOR M.  The  numerical  method  used  to  compute  the 
continuous  Fourier  transform  is  verified  by  applying  it  to  a  test  case 
which  is  solved  analytically  in  the  class  notes  (Ref  7:122).  The 
electric  field  is  written  as: 

E(t)  =  EQ*[EXP(-at)-EXP(-6t)]  (A.1) 

where:  =1x10^  volts/raeter 

a  =1x10^  second"^ 

6  =5x10^  second"^ 


This  model  describes  a  pulse  which  rises  to  a  peak  in  a  few  shakes  and 
then  decays  in  several  shakes  (Figure  A.I).  The  analytic  solution  is 
written  as: 


2’!  ( )  ( 3^+“'^ ) 


where:  - (  j)  -energy  ( joule* second/neter^) 

w  -angular  frequency  (radians/second) 


^  =  ^'o^  ^o 


The  method  used  to  perform  the  Fourier  transform  is  described  in  Chapter 
IV.  The  electric  field  is  truncated  at  shakes  and  the  Fourier 
transform  is  performed  on  the  resulting  pulse.  The  truncated  form  of 
the  pulse  is  shown  in  Figure  A.I. 


C  FIELD 


The  solution  from  EMPFRE  is  compared  to  the  analytic  solution  in 
Figure  A. 2.  At  low  frequencies,  EMPFRE  underestimates  the  energy 
because  of  the  pulse  truncation.  At  frequencies  above  100  MHz,  EMPFRE 
overestimates  the  energy  because  of  the  sampling  frequency.  The 
parametric  studies  deal  with  the  frequency  range  between  100  KHz  and  100 
MHz.  The  sampling  rate  was  chosen  to  provide  accurate  results  up  to  100 
MHz.  Accuracy  beyond  100  MHz  would  require  higher  sampling  rates. 

ELECTRIC  FIELDS.  HEMPB  is  a  computer  code  used  at  the  Air  Force 
Weapons  Laboratory  to  model  EMP.  Lt  D.  Benzer  provided  HEMPB  output  for 
use  as  a  test  case.  The  burst  parameters  are  provided  below. 


HOB  =  100  Km  (A. 3) 

Gamma  Yield  =  0.01  Ktons  (A. 4) 

Target  Location:  LONGITUDE  5°W 
LATITUDE  50°N 

Burst  Location;  West  of  Target 

Slant  Range  s  628  Km  (A. 5) 

Geomagnetic  field  r  5.6x10“^  Weber  s/m^  (A. 6) 

Gamma  Ray  Energy  =  2.5  MeV  (A. 7) 

Alpha  =  5.0x10®  second"^  (A. 8) 

Beta  =  0.55x10®  second”^  (A. 9) 

TP  s  3«0x10“®  seconds  (A. 10) 


The  peak  electric  and  magnetic  fields  are  computed  for  26  ranges  along 
the  slant  range  (Table  A.1).  Computations  start  at  a  range  of  257  Km 
from  the  burst  and  step  toward  the  target.  Computations  end  when  the 
Compton  current's  effect  on  the  EMP  become  small.  The  electric  field  is 
then  calculated  at  the  target  allowing  for  1/r  attenuation. 


JLE>^5EC0ND/I1ETER^>^2) 


The  burst  parameters  for  EMPFRE  were  generated  to  match  the  test 
scenario.  Figures  3.2  and  3.3  are  used  to  find  the  declination  and  dip 
angle  of  the  geomagnetic  field. 


Declination  =  15°  West  (A. 11) 
Dip  Angle  s  70°  North  'A. 12) 
X  =  (628  Km)  •  C0S(15°)  =  6.06E5  meters  (A. 13) 
Y  =  (628  Km)  •  SIN(15°)  =  1.63E5  meters  (A.I^) 
Z  =  1.0x10^  meters  <'A.15) 
Beta  =  5.55x10®  seconds"^  (A.16) 


Other  parameters  follow  directly  from  above. 

The  peak  electric  field  predicted  by  EMPFRE  is  285  vo  1  ts/'meter. 
The  prediction  from  HEMPB  is  198  volts/meter.  The  most  likely  cause  for 
the  error  is  due  to  the  flat  earth  approximation  and  its  effect  on  the 
atmosphere  profile  along  the  slant  range.  This  source  of  error  will 
have  no  impact  on  the  results  of  the  parametric  studies  because  the 
targets  are  at  ground  zero. 
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ALT(Km) 

RANGE(Km) 

E-MAX(V/M) 

1 

60.000 

256.80 

0 

2 

57.910 

272.26 

70.744 

3 

55.857 

287.73 

102.61 

4 

53.841 

303.19 

135.70 

5 

51.851 

318.66 

170.13 

6 

49.917 

334.13 

209.13 

7 

48.011 

3^9.59 

[ 

250.58 

8 

46.142 

i  365.06 

i  294.02 

1 

9  i 

44.309 

:  380.53 

!  1 

338.72 

10  i 
1 

42.512 

1  395.99  ' 

333.57 

1 

i 

40.752 

1 

1  411.46 

1  419.32 

12 

i  39.030 

1  426.92 

1 

j  444.98 

13 

37.344 

j  442.39 

j  460.45 

1 

14 

35.695 

457.86 

i  467.27 

15 

'  34.083 

i  473.32 

i  464.87 

16 

32.508 

488.79 

455.48 

17 

1  30.970 

i  504.25 

442.72 

18 

i 

1  29.469 

! 

519.72 

!  429.48 

19 

28.005 

i 

1  535.19 

i 

20 

26.578 

550.65 

i 

21 

25.188 

j  566.12 

1 

1 

i 

0 

1 

627.98 

.  198.45 

Table  A.1  Maximum  Electric  Fields  From  HEMPB 
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APPENDIX  B 

PROGRAM  DESCRIPTION  AND  LISTING 

EMPFRE  calculates  the  gararaa  ray  induced  EMP  from  a  high  altitude 
nuclear  detonation.  The  program  is  written  in  FORTRAN  IV.  A  listing  of 
EMPFRE  is  provided  along  with  comments  to  identify  key  variables  and 
procedures.  A  detailed  explanation  of  the  methods  employed  is  provided 
in  Chapter  III.  Parameter  variations  are  made  by  editing  the  code.  The 
results  are  written  to  a  tape.  Another  program  accesses  the  tapes  and 
generates  a  plot  using  DISSPLA8.  EMPFRE  generates  the  electric  field 
and  continuous  Fourier  transform  for  one  scinario.  Execution  time  on 
the  CyBEfi-750  is  roughly  5  seconds  per  electric  field  evaluation  made. 

Key  flow  diagrams  are  given  in  "igures  3.i  and  B.T. 
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INITIALIZE  KEY  PARAMETERS  ! 
t 

1  SUBROUTINE  DYNAMIC;  CALCULATES 

I  COMPTON  ELECTRON  PARAMETERS 

f! 

L. 

y 

SETS  LIMITS  FOR  SPATIAL  INTEGRATION 

u 

NORMALIZE  YIELD  FUNCTION  , 

1 

i 

SUBROUTINE  EFIELD;  CALCULATES 

1  : 

1 

[ 

ELECTRIC  FIELDS 

( 

SUBROUTINE  FOURIER:  CALCULATES 
CONTINUOUS  FOURIER  TRANSFORM 


/  RECORDS  RESULTS  ON  TAPE  / 
C - - 

■_EN^/ 


FIGURE  B.1  Flow  Diagran  For  Program  EMPFRE 
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FIGURE  B.2  Flow  Diagram  For  Subroutine  EFIELD 
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